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Abstract. We study the problem of two interacting particles in a one-dimensional quasiperiodic lattice of 
the Harper model. We show that a short or long range interaction between particles leads to emergence of 
delocalized pairs in the non-interacting localized phase. The properties of these Freed by Interaction Kinetic 
States (FIKS) are analyzed numerically including the advanced Arnoldi method. We find that the number 
of sites populated by FIKS pairs grows algebraically with the system size with the maximal exponent 
b = 1, up to a largest lattice size N = 10946 reached in our numerical simulations, thus corresponding to 
a complete delocalization of pairs. For delocalized FIKS pairs the spectral properties of such quasiperiodic 
operators represent a deep mathematical problem. We argue that FIKS pairs can be detected in the 
framework of recent cold atom experiments [M. Schreiber et al. Science 349 , 842 (2015)] by a simple setup 
modification. We also discuss possible implications of FIKS pairs for electron transport in the regime of 
charge-density wave and high T, superconductivity. 

PACS. 05.45.Mt Quantum chaos; semiclassical methods - 72.15.Rn Localization effects (Anderson or 
weak localization) - 67.85.-d Ultracold gases 


1 Introduction 

The Harper model ^ describes the quantum evolution of 
an electron in a two-dimensional periodic potential in a 
magnetic field. Due to periodicity it can be reduced to 
a one-dimensional Schrodinger equation on a quasiperi¬ 
odic lattice known as the almost Mathieu operator. This 
equation is characterized by a dimensional Planck con¬ 
stant determined by the magnetic flux through the lat¬ 
tice cell. The complex structure of the spectrum of this 
model was discussed in [2] and was directly demonstrated 
in [3]. As shown by Aubry and Andre [4], for irrational 
flux values a/27r this one-dimensional (ID) system has a 
metal-insulator transition with ballistic states for A < 2 
(large hopping) and localized states for A > 2 (small hop¬ 
ping). The rigorous proof is given in [^. The review on 
this model can be found in [6] and more recent results are 
reported in wm- 

It is interesting to study the case of Two Interact¬ 
ing Particles (TIP) in the Harper model. The model with 
Hubbard interaction between two particles was introduced 
in and it was shown that an interaction of moderate 
strength leads to the appearance of a localized component 
in the metallic non-interacting phase at A < 2 while in the 
localized phase A > 2 such an interaction does not signif¬ 
icantly affect the properties of localized states. Further 
studies also showed that the interactions provide only an 
enhancement of localization properties mm- 

These results for the Harper model show an opposite 
tendency compared to the case of TIP in the ID Anderson 
model with disorder where moderate Hubbard interaction 


leads to an increase of the localization length for TIP com¬ 
paring to the non-interacting case [iiiiiiiiiiiiiiiinj. 

Thus the result of Flach, Ivanchenko, Khomeriki [TS] 
on appearance of delocalized TIP states at certain large 
interactions in the localized phase of the Harper model 
at A > 2 is surprising and very interesting. In a certain 
way one has in this TIP Harper model the appearance of 
Freed by Interaction Kinetic States (FIKS). In this work 
we investigate the properties of these FIKS pairs in more 
detail using numerical simulations for the time evolution 
of wave functions and a new approach which allows to de¬ 
termine accurate eigenvectors for large system sizes up to 
~ 10"*^ (corresponding to a two-particle Hilbert space of 
dimension ~ 10®). This approach is based on a combina¬ 
tion of the Arnoldi method with a new, highly efficient, 
algorithm for Green’s function evaluations. 

We note that the delocalization transition in the Harper 
model has been realized recently in experiments with non¬ 
interacting cold atoms in optical lattices |19j . Experiments 
with interacting atoms have been reported in |20| and 
more recently in [21] showing delocalization features of 
interactions. Thus the investigations of the properties of 
FIKS pairs are of actual interest due to the recent ex¬ 
perimental progress with cold atoms. We will discuss the 
possible implications of FIKS pairs to cold atom and solid 
state experiments after presentation of our results. 

The paper is composed as follows: we describe the 
model in Section]^ the new Green function Arnoldi method 
is introduced in Section the analysis of time evolution 
of wave functions is presented in Section [^ the proper- 
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ties of FIKS eigenstates for the Hubbard interaction are 
described in Section and for the long rang interactions 
in Section properties of FIKS eigenstates in momen¬ 
tum and energy representations are analyzed in Section!? 
possible implications for the cold atom experiments [20 
ini] are discussed in Section]^ the dependence on the flux 
parameter is studied in Section [^ and the discussion of the 
results is presented in Section [10| 

2 Model description 

We consider particles in a one-dimensional lattice of size 
N. The one-particle Hamiltonian for particle j is given 
by: 

= r(^) (1) 

tO) =<x + l\j + h. c^, (2) 

X 

V^^'> =^Vi{x)\x>j<x\j. (3) 

X 

The kinetic energy is given by the standard tight- 
binding model in one dimension with hopping elements 
t = —1 linking nearest neighbor sites with periodic bound¬ 
ary conditions. We consider a quasiperiodic potential of 
the form Vi{x) = Acos(aa: -I- /3) which leads for A > 2 
to localized eigenfunctions with localization length i = 
l/log(A/2) [3]. Usually one chooses a = 27r(-\/5 — l)/2 
such that a/{2'K) « 0.61803 is the golden ratio, the “most” 
irrational number. For time evolution we manly use the 
golden mean value (together with the choice /3 = 0) while 
for the eigestates we mainly use the rational Fibonacci ap- 
proximant a —> 27r/„_i//„ where /„ is a certain Fibonacci 
number and where the system size is just N = fn- Further¬ 
more, in order to avoid the parity symmetry with respect 
to X— >-N — X at 13 = 0 (that leads to an artificial eigen¬ 
value degeneracy) we choose for this case (3 = — l)/2. 

We will see later that this Fibonacci approximant of a 
is very natural and useful in the interpretation at finite 
system sizes (especially with respect to Fourier transfor¬ 
mation). In our main numerical studies for the eigenvec¬ 
tors we consider system sizes/Fibonacci numbers in the 
range 55 < /„ < 10946 and the parameter A is always 
fixed at A = 2.5 with a one-particle localization length 
£= l/log(A/2) « 4.48 jHIH] . 

In Secs. [8] and i we also consider different irrational 
values of a/(27r) (or suitable rational approximants for 
finite system size). This is motivated by the recent exper¬ 
iments of Ref. |5T] and interest to the overall dependence 
of the FIKS properties on the flux parameter a. 

We now consider the TIP case, when each particle is 
described by the one-particle Hamiltonian h^3) ^ is cou¬ 
pled by an interaction potential U(xi — X 2 ) with another 
particle. Here we use U{x) = C//(l -|- w|a:|) for |a:| < Ur 
[22] and U{x) = 0 if \x\ > Ur with Ur being the interac¬ 
tion range, U is the global interaction strength and w is 
a parameter describing the decay of the interaction. We 
choose mostly w = 0 but in certain cases also rc = 1. The 


case Ur = 1 corresponds to the case of the on-site Hub¬ 
bard interaction studied in msi- Here we consider both 
symmetric two-particle states (bosons) and (for Ur > 2) 
also anti-symmetric two-particle states (fermions). 

The total two-particle Hamiltonian is given by 

H = + /i(2) + u (4) 

where 

U = ^ U{xi - X2)\xi,X2><Xi,X2\ (5) 

Xl ,X2 

is the interaction operator in the two-particle Hilbert space 
and with the notation \xi,X 2 >= \xi >1 \x 2 >2 for the 
non-symmetrized two-particle states. 

Our aim is to determine if the interaction may in¬ 
duce at least partial delocalization, i. e. at least for some 
eigenstates at certain energies. This can be done by a 
time evolution calculation from the Schrodinger equation 
using a Trotter formula approximation (see Sec. or 
by a numerical computation of (some) eigenfunctions of 
H. The size of the (anti-)symmetrized Hilbert space is 
N 2 = N{N + s)/2 « N‘^j2 with s = 1 (s = —1) for the 
boson (fermion) case and therefore a direct full numerical 
diagonalization of H is limited to N smaller than a few 
hundred, e.g. N < 250 [18] . 

Since the Hamiltonian H corresponds to a sparse ma¬ 
trix one can in principle apply the Arnold! method [231 
1241125] or more precisely, since H is a Hermitian matrix, 
the Lanczos method [26] . to determine certain eigenval¬ 
ues and eigenvectors. In the next Section, we will present 
a new method based on the particular structure of H, 
the Green function Arnoldi method, which is even more 
efficient than the standard implicitely restarted Arnoldi 
method. Thus, it allows to study larger system sizes, to 
obtain more eigenvalues, for much more parameter val¬ 
ues and with virtually exact eigenvalues and eigenvectors, 
i. e. d'^E{iP) ~ 10 - 28 - 10 - 2 O implying that there are only 
numerical rounding errors due to the limited precision of 
standard double precision numbers. The description of the 
Arnoldi method and definition of 6‘^E(ip) are given in Ap¬ 
pendix [^ 

3 Green's function Arnoldi method 

Let E be some energy value for which we want to deter¬ 
mine numerically eigenvalues of El close to E and the cor¬ 
responding eigenvectors. Furthermore let G = {E — H)~^ 
be the Green function or resolvent of H at energy E. The 
idea of the Green function Arnoldi method is to apply the 
Arnoldi method to the resolvent G and not to H which 
is sufficient since the eigenvectors of G are identical to 
those of H and the eigenvalues Ej of H can be obtained 
from the eigenvalues 7 ^ of G simply by Ej = E — l/ 7 j. 
The important point is that the largest eigenvalues jj of 
G, which result from the simple Arnoldi method, provide 
exactly the eigenvalues Ej close to a given value E which 
we may choose arbitrarily. Therefore it is not necessary 
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to apply the quite complicated (and rather expensive) im¬ 
plicitly restarted Arnold! method in order to focus on a 
given energy interval. 

For this we need an efficient method to evaluate the 
product G\ip > of G to an arbitrary vector \ip > and an 
arbitrary value of E. We have developped a new, highly 
efficient, numerical algorithm to determine G\lp> with a 
complexity 0{U^N^) for an initial preparation step at a 
given value of E and 0{N^) for the matrix vector multi¬ 
plication, provided the value of E is kept fixed. For larger 
system sizes, when localization of one-particle eigenstates 
can be better exploited, the complexity of the matrix vec¬ 
tor multiplication can even be reduced to 0{cN‘^) with 
c ~ 10^ being a rather large constant. For comparison 
we remind that a naive matrix vector multiplication has 
a complexity of 0(7V|) = C>(7V^) assuming that the full 
matrix G has been calculated and stored previously. 

Our algorithm is based on the following “magic” exact 
formula: 

G = Go + Goil-UGo)-^UGo ( 6 ) 

where Gg is the resolvent at vanishing interaction and 
Go is its projection on the smaller subspace of dimension 
« UrN of sites in two-particle space where the interac¬ 
tion operator has a non-vanishing action. The computa¬ 
tion of Go and the matrix inverse in § can therefore be 
done with 0{U^N^) operations and has to be done only 
once for a given value of the Green function energy E. 
The full matrix Go does not need to be computed since 
we can efficiently compute the product Go\^> on a given 
vector I if > using a transformation of | (/? > from position 
to energy representation (in the basis of non-interacting 
two-particle product eigenstates) where Gq is diagonal 
and a further transformation back to position represen¬ 
tation. Both transformations can be done with complex¬ 
ity 0{N^) due to the product property of non-interacting 
two-particle eigenstates. Therefore allows to compute 
the product G\tf> also for the full resolvent G with 0{N^) 
operations which is exactly what we need to apply the 
Arnold! method to G. A second, even more efficient, vari¬ 
ant of the Green function Arnold! method actually uses 
directly vectors in energy representation thus reducing the 
number of necessary transformation steps by a factor of 
two and also provides certain other advantages. These and 
other details of this approach are described in Appendix 
while Appendix provides the proof of (|^. 


vector \ip{t) >, then transforming the resulting vector to 
momentum representation by Fast Fourier Transform us¬ 
ing the library FFTW [ 30 ], applying exp{—iHpAt) (di¬ 
agonal in momentum representation) and finally retrans¬ 
forming the vector back to position representation. For a 
finite value of At Q can be viewed as the “exact” time 
evolution of a “modified” Hamiltonian with E[ corrected 
by a sum of (higher order) commutators of E[p and Et^. 
We have chosen At = 0.1 and verified that it provides 
quantitatively correct results for the delocalization prop¬ 
erties and its parameter dependence (this was done by 
comparison with data at smaller At values). This integra¬ 
tion method for the time evolution already demonstrated 
its efficiency for TIP in a disordered potential [T3] . 

In all our numerical studies we fix A = 2.5 which has 
a modest one-particle localization length mM- The main 
part of studies is done for the irrational golden value of 
flux or rotati on n umber a/{2TT) = {'/b— l)/2 (all Sections 
except Secs. For the time evolution we choose the 

quasimomentum at /3 = 0 and use the system size N = 512 
with an initial state with both particles localized at the 
center point xo = N/2 with |'!/'(0) >= |a:o,a:o > for the 
boson case or an anti-symmetrized state with one-particle 
at position xo and the other one at position xq — 1, i. e. 
|■^/'(0) >= (|a;o,(a;o - I) > -|(a:o - I),a;o >)/V^, for the 
fermion case. 

To study the localization properties we use the one- 
particle density of states: 


Pi{x) = '^\ <X,X2 |V’> 


( 8 ) 


representing the probability of finding one-particle at po¬ 
sition X. We are interested in the case where only a small 
weight of density is delocalized from the initial state. Thus, 
we introduce an effective one-particle density without the 
20 % center box by using Pes{x) = Gpi{x) for 0 < a; < 
0.4iV or Q.QN < x < N and Pes{x) = 0 for 0.4A^ < x < 
0.6iV. Here G is a constant that assures the proper nor¬ 
malization Pes(x) = 1. Using this effective density we 
define two length scales to characterize the (low weight) 
delocalization which are the inverse participation ratio 


6pr — 



-1 


(9) 


4 Time evolution 

We start our numerical study with a calculation for the 
time evolution with respect to the Hamiltonian Q using 
a Trotter formula approximation: 

\'ip{t -1- At) >= exp{—iEIpAt) exp{—iE[xAt) \'4’{t) > (7) 

with Hp = -k T(2) and -k + jj_ The 

time evolution step Q is valid for the limit of small At 
and allows for an emcient evaluation by first applying 
exp{—iE[xAt) (diagonal in position representation) to the 


which gives the approximate number of sites over which 
the density (outside the 20% center box) extends and the 
variance length {(x — with 

{{x-xo)'^) =^{x-xo)‘^ Pes{x) . (10) 

X 

Fig.[l] shows the dependence of both length scales on 
the interaction strength U for values up to C/ < 20 and dif¬ 
ferent cases of interaction range Ur and decay parameter 
w at iteration time t = 5120 (or t = 20480 for the boson 
case with Ur = 7 and w = 1). For each case there are a few 
values of interaction strength where the delocalization is 
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Fig. 1. Inverse participation ratio ^ipr and variance length 
= {{x — xo)^)^^'^ of the time evolution two-particle 
state for system size N = 512 and iteration time t = 5120 (or 
t = 20480 for bottom left panel) versus interaction strength U. 
The initial state at t = 0 is localized either with both parti¬ 
cles in the center position xq = N/2 (boson case) or antisym¬ 
metrized with one-particle in position xq and the other par¬ 
ticle in position xo — 1 (fermion case). Both quantities have 
been calculated from an effective one-particle density with¬ 
out a center box of size 20% (with respect to system size). 
The different panels correspond to different cases of interac¬ 
tion range Ur, decay parameter w and boson/fermion case. 
Here a/(2Tv) = (\/5 — l)/2 and /3 = 0. 


rather strong, even if the weight of the delocalized compo¬ 
nent is relatively small. For the Hubbard interaction case 
Ufi we find the two interesting values U = 4.5 and U = 7.4 
in a rather good agreement with the results of Ref. [18]. 
However, a closer inspection of the one-particle density re¬ 
veals that there is still a strong localized main peak close 
to initial point Xq and the delocalization only applies to 
a small weight of the initial state. We also note that the 
quantity ( ^captures peaks in {7 in a more clear way com¬ 
pared to (10). We attribute this to additional fluctuations 
added by aHarge distance from Xq to x values outside of 
the central box. 

The localized main peak can be understood by the as¬ 
sumption that only a small fraction of (two-particle) eigen¬ 
vectors with specific energy eigenvalues are delocalized 
while the other eigenvectors remain strongly localized. In¬ 
deed, the initial vector |V'(0) >, localized at xq and ex¬ 
panded in a basis of two-particle energy eigenstates, con¬ 
tains contributions from all possible energy eigenvalues. 
The time evolution from the Schrodinger equation only 
modifies the phases of the energy expansions coefficients 
but not the amplitudes and therefore the wave packet at 
arbitrary time |'!/'(^) > contains rather uniform contribu¬ 
tions from the same energy values. Obviously the delo¬ 
calization effect in the wave packet only happens for the 
small weight corresponding to the limited fraction of de¬ 
localized eigenvectors while the other contributions form 
the central peak close to the initial position. 


Table 1. Time evolution parameters for certain cases of short 
and long range interactions for interaction values with strong 
delocalization. All rows except the last one correspond to the 
boson case and the last row to the fermion case. The iteration 
time is t = 5120 except for the case with U = 16.9, Ur = 7 
and w = 1 where t = 20480. Here a/(27r) = (\/5 — l)/2 and 
/3 = 0 . 


u 

Ur 

w 

?IPR 

(H) 



4.4 

1 

0 

129.16 

-3.0756 

0.2257 

0.04175 

4.5 

1 

0 

125.22 

-3.0645 

0.2454 

0.0383 

4.7 

1 

0 

148.56 

-3.0347 

0.2594 

0.02596 

7.2 

1 

0 

109.53 

1.8072 

0.4891 

0.05801 

7.4 

1 

0 

136.60 

1.1369 

3.0897 

0.04102 

7.8 

1 

0 

15.13 

1.8151 

0.6851 

0.0001974 

8.0 

5 

0 

89.26 

8.7256 

0.3260 

0.01406 

16.9 

7 

1 

136.06 

10.1893 

0.5026 

0.03268 

10.9 

5 

0 

243.17 

10.8879 

0.4431 

0.0795 


We have therefore computed a tail state liptaiiit) > from 
the wave packet \ip(t) > by removing (putting to zero) a 
big 60% center box in a similar way as for PeB{x) (but 
in the two-particle space and using a larger center box). 
The energy eigenvectors who contribute to |V’taii(i) > ob¬ 
viously only cover the delocalized eigenvectors and assum¬ 
ing that the latter exist only for certain specific energies 
we can try to determine this energy range (for delocaliza¬ 
tion) by computing the expectation value (H) of H and 
its energy variance [see Eq. ( j^ ] with respect to jV'taiKO > 
(after proper renormalization of |'0taii(i) >)■ Furthermore 
the square norm jjV'taiilOlP) which is the probability of 
propagating outside the 60% centerbox, gives also a good 
measure for the delocalization effect. 

In Table we show for certain cases with strong de- 
localization the values of the quantities ^ipr, {H), 5‘^E 
and ||'0taii(i)|P- For Ur = 1 and the first peak at f7 « 4.5 
the maximum for ^ipr corresponds to U = 4.7 while the 
maximum of ||'0taii(i)|P corresponds toU = 4.4. Therefore 
the intermediate value U = 4.5 used in Ref. [18] is indeed 
promising. For all these three values of U the average en¬ 
ergy (H) « —3.05 of the tail state corresponds rather 
well to the approximate eigenvalue region E « —3.1 at 
U = 4.5 for delocalized eigenstates found in [TB] and con¬ 
firmed by our detailed eigenvector analysis presented in 
the next Section. Furthermore, the corresponding energy 
variance is indeed rather small. 

For Ur = 1 there is also a second local maximum of 
CiPR at 17 = 7.4 and close to this value there is also a local 
maximum of ||V'taii(l)||^ at 17 = 7.2. We have also included 
in Table [T] the value U = 7.8 which is close to the second 
interaction value U = 7.9 used in Ref. m- The value 
U = 7.8 seems less optimal but our eigenvector analysis 
shows that this value is quite optimal for two different en¬ 
ergy ranges E « 1.8 and E « —2.8 with well delocalized 
eigenstates for both energies. According to Table the 
average energy of the tail state is {H) « 1.8 for U = 7.2 
and U = 7.8 but with a somewhat larger value of the 
variance (in comparison to the case U = 4.5) indicating 
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that the main contributions in the tail state arise from the 
first energy range E Ki 1.8 but the second value E « —2.8 
provides also some smaller contributions therefore increas¬ 
ing the variance. For U = 7.4 the average energy of the 
tail state is even reduced to {H) « 1.1 and the variance 
S^E « 3.1 is quite large which indicates clearly that for 
this case both energy ranges have more comparable contri¬ 
butions in the tail state. In Fig. 1(a) of Ref. |TH] these two 
energy values can be roughly identified with a somewhat 
stronger delocalization at i? « —2.8. Our eigenvector cal¬ 
culations (see next Section) for larger system sizes confirm 
that for modest values of system sizes the delocalization is 
stronger at E « —2.8 but at larger sizes it is considerable 
stronger at E « 1.8. 

The values of HV'taiKOlP between 10“'* and 5.8 x 10“^ 
represent the weight of the delocalized eigenstates in the 
wave packets. These values are significantly smaller than 
unity showing that the main contribution still corresponds 
to the central peak at xq and the localized eigenstates at 
other energy values but they are also considerably larger 
than the values ^ 10“*'* for U values with minimal (or 
absent) small weight delocalization. In general, the maxi¬ 
mal values of U for the two length scales shown in Fig. 
correspond rather well also to the local maximal values for 
||V'taii(0lP- For the other three cases of Fig. with long 
range interaction we can also identify certain values of U 
with rather strong delocalization (for both length scales 
and the squared tail norm). According to Tablewe find 
for these three cases ^ipr ~ 10 ^, ||V'taii(0lP ~ 
rather sharp average energy values of the tail state with a 
small variance. 

We have repeated this type of analysis also for many 
other long range interaction cases and in certain cases we 
have been able to identify optimal values of U and E for 
strong delocalization where the approximate energy ob¬ 
tained from the time evolution tail state was used as ini¬ 
tial value of E for the Green function Arnoldi method to 
compute eigenstates (see Sec.[^. 

We also computed the inverse participation ratio and 
the variance length using the full one-particle density of 
states (including the center box) and also these quanti¬ 
ties have somewhat maximal values at the optimum U 
values for delocalization found above but their maximum 
values are much smaller than the length scales shown in 
Fig. 0 Therefore it would be more difficult (or impossible) 
to distinguish between small weight long range delocaliza¬ 
tion and high weight small or medium range delocalization 
(i.e. where the full wave packet delocalizes but for a much 
smaller length scale). For this reason we prefer to compute 
the inverse participation ratio and the variance length us¬ 
ing the effective one-particle density without center box 
and with the results shown in Fig. 

In Fig. we show the density plots of a zoomed re¬ 
gion of the time evolution state for the four cases of Fig. 
[h and the optimal delocalization values for U {U = 4.5 
for Ur = 1 and the three values given in Table for 
the cases with Ur > 1 and also mentioned in the fig¬ 
ure caption of Fig. [^. The zoomed region correspond to 
a box of size 205 x 205 with left bottom corner at posi- 



Fig. 2. Density plot of time evolution state for t = 5120 (or 
t = 20480 for bottom left panel), system size N = 512, the four 
cases of Fig. [^and with a value of U corresponding to strongest 
delocalization: U = 4.5, Ur — 1, boson case (top left panel), 
U = 8, Ur = 5, w = 0, boson case (top right panel), U = 16.9, 
Ur = 7, w = 1, boson case (bottom left panel), U = 10.9, 
Ur = 5, w = 0, fermion case (bottom right panel). We show 
only zoomed region of size 205 x 205 with left bottom corner 
at position xi = X 2 = 307 which corresponds to the right/top 
boundary of the 20% center box. The colors indicate red for 
maximum, green for medium and blue for minimum values 
(same distribution of colors in other figures of density plots). 


tion xi = X 2 = 307. This value corresponds exactly to 
the right/top boundary of the 20% center box which has 
been removed when determining the effective one-particle 
density of states Peg{x). For positions inside the center 
box between 205 and 306 the time evolution state has a 
strong peaked structure with considerably larger values 
of the amplitude than the right/top part shown in Fig. 

The left/lower part (between 0 and 204) is similar in 
structure with similar amplitudes to the right/top part. 
Fig. [^clearly confirms the complete small weight delocal¬ 
ization along the diagonal xi « X 2 of the wave packet at 
sufficiently long iterations times t = 5120 (or t = 20480 
for the case with Ur = 7 and w = 1). 

The time evolution of the one-particle density of states 
can be seen in Fig. with its time dependence correspond¬ 
ing to the vertical axis and position dependence corre¬ 
sponding to the horizontal axis for the same cases and 
parameters of Fig.[^ In all cases one can identify a strong 
central peak at xq and a low weight delocalization with a 
characteristic length scale increasing linearly in time, thus 
corresponding to a ballistic dynamics already observed for 
the Hubbard interaction case in [TB] . One can also observe 
in Figs. and 1^ that for U = 8.0, Ur = 5, w = 0, boson 
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Fig. 3. Density plot for the time dependence of one-particle 
density from the time evolution state with a:-position (0 < a; < 
512) corresponding to the horizontal axis and time t {0 < t < 
5120 or 0 < t < 20480 for bottom left panel) corresponding 
to the vertical axis. The four panels correspond to the same 
parameter values of U, Ur, w and boson/fermion cases as in 
four panels of Fig. 


case, the weight of the delocalized part of the wave packet 
is minimal of the four shown cases which is in agreement 
with the lowest value of ||^/’taii(t)|P for the same case. 


5 Eigenstates for Hubbard interaction 


In this Section we present our results for the two-particle 
eigenstates for the case of the Hubbard interaction with 
Ur = 1. In order to characterize the delocalization proper¬ 
ties of eigenstates we use two quantities. One is the inverse 
participation ratio in position representation obtained 
from the one-particle density of states ^ of eigenstate 
IV’>, by 



Another one is the inverse participation ratio in energy 
representation obtained from an expansion of a two- 
particle eigenstate \tjj> ot H in the basis of non-interacting 
energy product eigenstates (of Hq) by 


= 






-1 


( 12 ) 


The quantity is identical to the “participation num¬ 
ber” used in Ref. [I8]. It is similar (but different) to the 


quantity used in the previous Section, but for the full 
one-particle density and not the effective density without 
the 20% center box. Thus counts the number of x- 
positions over which the one-particle density extends and 
obeys the exact inequality ^x < N. It is not to be con¬ 
fused with the inverse participation ratio in the two par¬ 
ticle (xi,a; 2 )-space, a quantity we did not study. Instead 
we use the other quantity that counts the number of 
non-interacting energy product eigenstates of Hq which 
contribute in the eigenstate. This quantity may be larger 
than N as we will see for the case of long range interac¬ 
tions in the next Section. It is very convenient to deter¬ 
mine ^E with the second variant of the Green function 
Arnoldi method where the main computations are done 
in the energy representation using the non-interacting en¬ 
ergy product eigenstates as basis states. For the 

case of two particles localized far away from each other, 
the quantity ^e is very close to unity while ^x is closer to 
3-4 due to the finite localization length of the one-particle 
Harper problem. For a ballistic delocalized state along the 
diagonal xi = X 2 we expect that both ^x and ^e are ~ CN 
with some constant C of order or a bit smaller than unity. 

In this and the next Sections we choose the system 
size to be a Fibonacci number iV = /„, the rational case 
a/(27r) = fn-i/fn and /3 = (v^- l)/2. However, we have 
verified that the strong delocalization of eigenstates for 
certain values of U and E is also valid for the irrational 
case for arbitrary N with a/(27r) = {y/5 — l)/2 and /3 = 
0. For example for Ur = !,[/ = 4.5, E « —3.1 {U = 
7.8, E « —2.8) we find for the rational case with N = 
4181 that the eigenstate with maximal ^e corresponds to 
E = -3.09901, ^E = 795.960 and = 1172.887 {E = 
-2.78600, ^E = 501.321 and = 475.573) while for the 
irrational case with N = 4000 we have E = —3.09963, 
^E = 763.440 and ^x = 889.854 {E = -2.78716, ^e = 
559.130 and = 588.186). 

We consider as system size N all Fibonacci numbers 
between 55 and 10946. For each system size we apply the 
Green function Arnoldi method with a typical Arnoldi 
dimension ha ~ 0.7N-0.8N slightly smaller than N ex¬ 
cept for the largest case N = 10946 for which we choose 
UA = 2000 or HA = 3000 and the smallest cases TV = 55 or 
TV = 89 where we choose ua ^ 300-400. From all ua Ritz 
eigenvalues we retain only those with a minimal quality 
requirement of S‘^E('ip) < 10“® which corresponds roughly 
to 2/3 of all UA eigenvalues. It turns out that among these 
“acceptable” eigenvalues most of them are virtually exact 
with S'^E{tfj) < 10“^° (or even better), especially for the 
eigenvalues closest to the Green function energy E or with 
rather large values of ^e or ^x- Only some eigenvalues at 
the boundaries E ± AE (with AE depending on TV and 
Ua) of the obtained energy band were of modest quality 
with 6'^E{'tjj) between 10“^° and 10“®. 

Goncerning the interaction strength U and the approx¬ 
imate energy range E we present here the detailed results 
for the eigenvectors of four cases which are U = 4.5 com¬ 
bined with E = —3.1, U = 7.2 combined with E = 1.8 and 
also the less optimal interaction strength U = 7.8 with two 
possible energy values E = —2.8 and E = 1.8. For three of 
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theses cases {U = 4.5, U = 7.2 and U = 7.8 with E = 1.8) 
the approximate energy range can be obtained as the av¬ 
erage energy {H) of the tail state computed from the time 
evolution and given in Table [l] For the last case the sec¬ 
ond interesting energy value E = —2.8 for U = 7.8 can 
be found by exact diagonalization for small system sizes 
(iV = 55 and N = 89) and was also identified in Fig. 1(a) 
of Ref. [IH] . (Actually, the Green function Arnoldi method 
is for small system sizes also suitable for a full matrix diag¬ 
onalization by choosing ua = N{N -\- 1)/2 identical to the 
dimension of the symmetrized two-particle Hilbert space.) 

The Green function Arnoldi method requires to fix a 
preferential energy for the Green function which deter¬ 
mines the approximate energy range of computed eigen¬ 
values and eigenvectors. For this we use a refinement pro¬ 
cedure where at each system size N this energy is either 
chosen as the eigenvalue of the eigenstate with maximum 

obtained from the last smaller system size or, for the 
smallest system size N = 55, as one of the above given 
approximate energy values essentially obtained as the av¬ 
erage energy of the time evolution tail state. This system¬ 
atic refinement is indeed necessary if one does not want to 
miss the strongest delocalized states since the typical en¬ 
ergy width of “good” eigenvalues provided by the method 
decreases rather strongly with increasing system size, e. g. 
AE ~ 10-3 for N = 10946. 

In this way we obtained indeed the strongest delocal¬ 
ized states up to the largest considered system size. How¬ 
ever, for N = 10946 we added one or two additional runs 
at some suitable neighbor values for E which allowed us to 
obtain a more complete set of delocalized states. We also 
made an additional verification that overlapping states, 
obtained by two different runs at different E values, were 
indeed identical for both runs and did not depend on the 
precise value of E used in the Green function Arnoldi 
method provided that the eigenvalue of the overlapping 
eigenstate was sufhciently close to both E values. In gen¬ 
eral, if one is interested in an eigenstate which by accident 
is close to the boundary of the good energy interval and 
is therefore of limited quality, one can easily improve its 
quality by starting a new run with a Green function energy 
closer to the eigenvalue of this state. 

In Fig. 1^ we show density plots for the strongest de¬ 
localized eigenstates (in ^e) for the two cases U = 4.5, 
E « —3.1 and U = 7.8, E « —2.8 and the three smallest 
system sizes N = 55, N = 89 and N = 144. In all cases 
the eigenstate extends to the full diagonal along xi « X 2 
with a width of about 7 sites {U = 4.5) or about 15 sites 
{U = 7.8) with a quasiperiodic structure of holes or strong 
peaks. One can also identify some additional peaks with 
\xi — X 2 \ 20-30 which can be interpreted as a resonant 

coupling of the main state with some product state of 
non-interacting one-particle eigenstates with both parti¬ 
cles localized at some modest distance a bit larger than 
the one-particle localization length i « 4.48 and where 
the eigenvalue of the main state is very close to the total 
energy of the product state. 

In Fig. and Fig. the strongest delocalized states 
for N = 1597 {N = 10946) and the same values of U and 



Fig. 4. Density plot of FIKS eigenstates with maximal value 
for system size N = 55 (top panels), A = 89 (center panels), 
N = 144 (bottom panels), Ur = 1 and interaction strength 
U = 4.5 (left column) or U = 7.8 (right column). The corre¬ 
sponding energy eigenvalues and values for both types of in¬ 
verse participation ratios are; Top left: E = —3.10334, = 

22.756, = 30.794. Top right: E = -2.75868, ^e = 33.274, 

G = 24.901. Center left: E = -3.09588, ^e = 50.742, = 

49.867. Center right: E = -2.78575, Cs = 35.139, G = 28.198. 
Bottom left: E = 3.09966, f,E = 61.373, = 63.353. Bottom 

right: E = -2.78596, ^e = 56.210, G = 47.958. 


approximate energy as in Fig. [^are shown as full states 
(only for N = 1597) and with three zoomed regions of 
size 100 X 100 at three different positions on the diagonal 
(for N = 1597 and N = 10946). Again the eigenstates 
extend to the full diagonal size with a certain width and 
one can identify a a quasiperiodic structure of holes and 
peaks and some resonant couplings to product states of 
non-interacting one-particle eigenstates. Higher quality gif 
files for the full eigenstate of these (and some other) cases 
are available for download at m- 

Figs. |4][^ also show that, apart from the common fea¬ 
tures, with increasing system size the eigenstates seem to 
become “thinner”, i. e. the weight of the hole parts seems 
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Fig. 5. Density plot of FIKS eigenstates with maximal value 
of for system size N = 1597, Ur — 1 (both columns) and 
interaction strength U = 4.5, energy eigenvalue E = —3.09644, 
= 616.638, = 716.050 (left column) or 17 = 7.8, E = 

-2.78777, is = 330.269, ix = 355.236 (right column). The 
first row corresponds to the full eigenstates and the other rows 
correspond to zoomed regions of size 100 x 100 with bottom left 
corner at position *1 = a ;2 = 0 (second row), x\ = X 2 = 700 
(third row) and xi = 2:2 = 1400 (fourth row). 



Fig. 6. Density plot of FIKS eigenstates with maximal value 
of £,e for system size N = 10946, Ur — 1 (both columns) and 
interaction strength U = 4.5, energy eigenvalue E = —3.09749, 
Cb = 2099.806, ^x = 3105.529 (left column) or U = 7.8, E = 
-2.78707, ^E = 952.498, ^x = 1147.965 (right column). All 
panels correspond to zoomed regions of size 100 x 100 with 
bottom left corner at position a;i = X 2 = 0 (first row), Xi — 
X 2 = 5000 (second row) and xi = X 2 = 10000 (third row). 


to increase and the strength of peaks seems to decrease, 
especially for the case U = 7.8 and approximate energy 
E = -2.8. 

Fig. shows a zoomed region of size 100 x 100 roughly 
in the middle of the diagonal for strongest delocalized 
eigenstates for N = 1597 and N = 10946 and the two 
cases U = 7.2 and U = 7.8, both with the approximate en¬ 
ergy E — 1.8. Globally one observes in Fig.j^the same fea¬ 
tures as in the Figs.[5][^for the previous two cases but with 
a detail structure on the diagonal which is significantly dif¬ 
ferent, i. e. quite large width and different pattern for the 
quasiperiodic peak-hole structure. One observes that the 
eigenstates for U = 7.2 are very compact while for U = 7.8 
they are a bit less compact, with more holes, but also with 
additional small satellite contributions from product pair- 
states at distance « 20 from the diagonal. These satellite 
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Fig. 7. Density plot of FIKS eigenstates with maximal value of 
for E ^ 1.8, Ur = 1,U = 7.2 (left column) orU = 7.8 (right 
column) and system size N = 1597 (top panels) or N = 10946 
(bottom panels). All panels correspond to zoomed regions of 
size 100 X 100 with bottom left corner at position xi = X 2 = 
700 (top panels) or x\ = X 2 = 5000 (bottom panels). The 
corresponding energy eigenvalues and values for both types of 
inverse participation ratios are: Top left: E = 1.79597, = 

638.916, = 506.113, Top right: E = 1.81744, = 475.972, 

= 359.239, Bottom left: E = 1.79652, ^e = 5694.610, = 

4834.890, Bottom right: E = 1.81741, ^e ~ 2086.088, = 

1843.227. The figures obtained for other zoomed regions (on 
the diagonal) for these states are very similar and these four 
eigenstates extend clearly to the full diagonal xi « X 2 and all 
values with 0 < a;i < A^. 


contributions are absent at U = 7.2. Apart from this the 
pattern for both cases in Fig. is rather similar, i.e. the 
FIKS eigenstates for E ^ 1.8, and U = 7.2 or U = 7.8 be¬ 
long to the same family but obviously the value U = 7.2 
is more optimal with a compacter structure, larger values 
of and This is also in agreement with the discussion 
of the time evolution states in the previous Section. It is 
interesting to note that even for the case U = 7.8 with a 
modest squared tail norm « 2 x 10“^ (instead of 5 x 10“^ 
for U = 7.2, see Tablethere are very clear FIKS eigen¬ 
states and even at two different energy regions. 

We have also calculated eigenstates up to system sizes 
N = 2584 for the additional case U = 7.2 and E « —2.8 in 
order to verify if the second energy value is also interesting 
for U = 7.2. Here one finds also some FIKS eigenstates but 
of reduced quality if compared to U = 7.8 and E « —2.8, 
i. e. smaller values of and and for larger system sizes 
the eigenstates do not extend to the full diagonal, i. e. 
about 20-40% of the diagonal is occupied for N = 2584. 
For this additional case we do not present any figures. 



4000 


N=10946 + 


N=10946 + 

N=6765 X 


N=6765 X 

N=4181 * 


N=4181 * + 

N=2584 □ + 


N=2584 □ + 

N=1597 X 


N=1597 ^ 

N=987 0 XX : 

ujf 2000 

N=987 0 * 


1000 
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Fig. 8. Inverse participation ratio of eigenstates versus eigen¬ 
value energy for the system sizes N = 987, 1597, 2584, 4181, 
6765, 10946 and Ur = 1. The left column of panels correspond 
to the inverse participation ratio ^e in energy representation 
and the right column to the inverse participation ratio in po¬ 
sition representation. First row of panels correspond toU = 4.5 
and the energy region E « —3.098, second (third) row of pan¬ 
els correspond to U = 7.8 and E « —2.787 {E « 1.817) and 
fourth row of panels correspond to 17 = 7.2 and E « 1.79. 


In Fig. 1^ both types of inverse participation ratios ^e 
and ^x of eigenstates are shown as a function of the energy 
eigenvalue for all four cases (corresponding to Figs. IB 
with energies in the interesting regions and for the six 
largest values of the system size between 987 and 10496. 
Both quantities increase considerably with system size and 
the overall shape of the cloud of points seems to be similar 
for each value of N but with a vertical scaling factor in¬ 
creasing with N. The figures for ^e and ^x are rather simi¬ 
lar with somewhat larger (maximum) values for ^x (except 
for U = 7.2 where the maximum value of ^e is larger). For 
U = 4.5 the energy region of delocalized states extends 
from E « —3.103 to A « —3.092 and for N = 10496 two 
supplementary runs with Green’s function energy values 
shifted to the left {E = —3.104) and right {E = —3.094) 
from the center {E = —3.0977) were necessary to obtain 
a complete cloud of data points. For U = 7.8 and ap- 
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Fig. 9. Largest inverse participation ratio (for given values of 
N and approximate energy) of FIKS eigenstates versus sys¬ 
tem size in a double logarithmic scale using all Fibonacci 
numbers between 55 and 10946. Top (bottom) left panel corre¬ 
sponds to ?7 = 4.5 {U = 7.2) and the energy region E ~ —3.1 
{E « 1.8). Top (bottom) right panel corresponds to U = 7.8 
and E « —2.8 {E « 1.8). The blue line with stars corresponds 
to the inverse participation ratio in position representation, 
the red line with crosses to the inverse participation ratio 
in energy representation and the green line to the power law fit 
= aN*‘ with fit results shown in Table The shown energy 
values in the panels refer to the eigenvector with maximal 
for the largest system size. Note that for given values of N 
and approximate energy the eigenstates with maximal ^x and 
maximal ^e may be different. 


proximate energy E = —2.8 the main region of delocal¬ 
ized eigenstates extends from E « —2.788 to ill « —2.786 
with a secondary small region at if « —2.789. For the 
secondary region and N = 10946 also an additional run 
with a shifted Green function energy was necessary. For 
U = 7.8 and approximate energy £1 = 1.8 the main re¬ 
gion of delocalized eigenstates extends from E « 1.8172 
to £1 « 1.8186 also with a secondary small region at 
E « 1.814 and for this secondary region and N = 10946 
also an additional run with a shifted Green function en¬ 
ergy was necessary. 

For U = 7.2 and approximate energy E = 1.8 the main 
region of delocalized eigenstates extends from E « 1.793 
to £1 « 1.797. For this particular case one observes the 
absence of eigenstates with very small values of « 1 
and «3-4. We have verified, by choosing different val¬ 
ues of the Arnold! dimension tia and the Green function 
energy, that the absence of such states is stable with re¬ 
spect to different parameters of the numerical method. 
Apparently in this energy region there are no strongly lo¬ 
calized product states (of one-particle energy eigenstates) 
with a modest distance between the two particles such 
that there would be some contribution of them in the ini¬ 
tial state used for the Arnold! method. There may still be 
other product states in this energy region but with the 
two particles localized further away such that the Arnoldi 
method cannot detect them. 


Table 2. Approximate energy E (for largest system size) and 
results of the power law fit ^x = aN^ for the same cases and 
data sets as in Fig. 


u 

Ur 

E 

a 

b 

4.5 

1 

-3.097 

0.940 

± 

0.137 

0.882 

± 

0.021 

7.2 

1 

1.797 

0.375 

± 

0.054 

1.003 

± 

0.021 

7.8 

1 

-2.787 

1.878 

± 

0.380 

0.698 

± 

0.029 

7.8 

1 

1.817 

0.559 

± 

0.073 

0.887 

± 

0.019 


The scenario of strongly delocalized eigenstates for cer¬ 
tain narrow energy bands found in Ref. |18| is clearly 
confirmed also for larger system sizes up to = 10946. 
However, the maximum values of and do not scale 
always linearly with N as can be seen in Fig. which 
shows the dependence of maximum values of and 
for all four cases (of interaction strength and approximate 
energy) as a function of the system size TV in a double 
logarithmic scale. Note that in Fig. the data points 
for maximum ^e (for given values of N, U and approx¬ 
imate energy) may correspond to other eigenstates than 
for the data points for maximum ^x-, i- e. the maximum 
values for the two quantities are obtained at two differ¬ 
ent eigenstates. For example for U = 4.5 and N = 6765 
the eigenstate with maximum ^e corresponds to E = 
—3.09771, ^E = 1861.131, ^x = 2538.299 while the eigen¬ 
state with maximum ^x corresponds to £1 = —3.09749, 
^E = 1406.560, ^x = 2573.484, a state which ranks on 
the 5th position in the list of states with maximum val¬ 
ues for ^E- However, despite such particular cases the ap¬ 
pearance of large values for ^e (strong delocalization in 
one-particle energy representation) or ^x (strong delocal¬ 
ization in position representation) are rather well corre¬ 
lated which is obvious since the transformation from en¬ 
ergy to position representation corresponds somehow to a 
“smoothing” on the length scale of the one-particle local¬ 
ization length i « 4.48. 

The results of the power law fit ^x = aN^ using the 
data sets of Fig. are shown in Table For U = 4.5 
01 U = 7.8 (both energy ranges) the fit values of the ex¬ 
ponent b, which are either close to 0.9 or 0.7, seem to 
indicate a kind of fractal structure of the eigenstates since 
even for the largest system sizes the corresponding eigen¬ 
states extend to the full length of the diagonal xi « X 2 - 
Therefore the reduction of ^x with respect to a linear be¬ 
havior in N is due to the internal structure (appearance 
of more “holes”). This is also in agreement with our above 
observation that delocalized eigenstates seem to become 
thinner for larger systems sizes and this effect is strongest 
for the case U = 7.8, E « —2.8 which also corresponds to 
the smallest value of the exponent b = 0.698 among the 
three cases. However, for U = 7.2 the exponent is rather 
precisely unity and no fractal or increasing hole structure 
(with increasing system size) is visible in the FIKS eigen¬ 
states (see also Fig. [^. 
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6 Eigenstates for long range interaction 


We now turn to the case of long range interactions with 
Ufi>l. We remind that we consider a model where the 
particles are coupled by the interaction potential U{xi — 
X 2 ) with U{x) = U/{l-\-w\x\) for \x\ < Ur [52] and U{x) = 
0 if I a; I > Ur. For the decay parameter w we mostly choose 
tc = 0 (i. e. “no decay”) or for the boson case also w = 1 
(decay ~ |a:i — X2\~^ provided that \xi — X2\ < Ur). 

We considered many different cases with 2 < Ur < 7 
and one case with Ur = 20 and performed for each case 
a time evolution analysis as described in Sec. to find 
good candidates of the interaction strength U for strong 
delocalization. Using the tail state analysis we also ob¬ 
tained suitable approximate energy values to start the 
Green function Arnoldi method for the smallest system 
size fV = 55 we considered. Then we refined the Green 
function energy for larger system sizes in the same way 
as described above. In many cases (but not always) this 
procedure leads to a nice data set of well delocalized two- 
particle eigenstates for a given narrow energy band. In 
certain cases the refinement procedure gets trapped at a 
“wrong” energy, i.e. which is promising for a particular 
small system size but where the localization saturates at 
some medium value for for larger system sizes or is sim¬ 
ply less optimal than some other energy. In these cases it 
might be useful to manually select a different eigenvalue 
obtained from the last smaller system (e. g. for iV = 55 or 
N = 89) to force the refinement of energies into a direction 
of stronger delocalized states. 

We mention that for the larger values oi Ur the com¬ 
putational cost [~ (NUr)^] and the memory requirement 

(NUr)"^] of the initial preparation part of the Green 
function Arnoldi method is considerably increased and 
therefore we have limited for these cases the maximal con¬ 
sidered system size to TV < 1597. 

In Fig. 10 we show the strongest delocalized state (in 
^e) for N = 610 and the case U = 14.0, Ur = 20, w = 0, 
boson case (top panels) and the three cases with Ur > 1 
already presented in Figs. [I][^ of Sec. (second to fourth 
row of panels). Goncerning the case Ur = 7, U = 16.9, 
w = 1, bosons (of Sec.|^, it turns out that for the eigen¬ 
state analysis the interaction strength U = 17.0 is some¬ 
what more optimal than the case oi U = 16.9. Therefore 
we show in Fig. 10 (and other figures in this Section) the 


case oi U = 17.0 instead oi U = 16.9. For each case the 
left column panel of Fig. shows the full state and the 
right column panel a zoomed region of size 100 x 100 with 
bottom left corner at position xi = X 2 = 200 for a better 
visibility. 

The energy eigenvalues of the three boson states in 
Fig. [1^ A = 14.00502, E = 8.79607 01 E = 10.22864 
(top three rows of panels) correspond quite well to the 
approximate energies obtained from the tail state anal¬ 
ysis of the time evolution wave packet for the same (or 
very similar) parameters: (H) = 14.00247, (H) = 8.72561 
or (i7) = 10.18926 (see also Table 0. However for the 
fermion case (fourth row of panels with U = 10.9, Ur = 5, 
ic = 0) the energy eigenvalue of the strongest delocalized 
state at fV = 610 is if = 11.53294 while the approxi- 



Fig. 10. Density plot of FIKS eigenstates for different cases 
of long range interaction Ur > 1 with maximal inverse partici¬ 
pation ratio in energy representation for N = 610. The left 
column corresponds to the full eigenstate and the right column 
to a zoomed region of size 100 x 100 with bottom left corner 
at position xi = X 2 = 200. First row: Ur = 20, w = 0, U = 14, 
boson case, energy eigenvalue E = 14.00502, = 263.410, 

G = 350.519. Second row: Ur = 5, w = 0, U = 8, boson case, 
E = 8.79607, Cs = 787.137, G = 397.779. Third row: Ur = 7, 
w = 1, U = 17, boson case, E = 10.22864, ^e = 635.918, 
G = 307.585. Fourth row: Ur = 5, w = 0, U = 10.9, fermion 
case, E = 11.53294, ^e = 535.618, G = 360.478. 
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mate energy obtained from the tail state analysis (H) = 
10.88786 is somewhat different. Here the refinement pro¬ 
cedure to optimize leads already at the first Green 
function Arnoldi calculation for N = 55 and ua = 400 
to an energy shift from 10.9 (as initial Green’s function 
energy) to 11.5 (as eigenvalue of the eigenstate with max¬ 
imum ^e)- However, optimizing for (instead of ^e) or 
fixing manually the value E = 10.9 for N = 144 results in 
a different set of strongly delocalized eigenstates close to 
the energy E = 10.84 with somewhat smaller values for 
^E but larger values for than the first set of delocalized 
eigenstates at E = 11.53. 

The eigenstates shown in Fig. [^have the same com¬ 
mon features as the eigenstates shown in Figs. |4][7] for the 
Hubbard short range interaction discussed previously such 
as extension to the full diagonal at xi « X 2 -, a certain 
width of ^ 10-20 sites, quasiperiodic structure of holes 
and peaks etc. but the detail pattern is specific for each 
case. For the very long interaction range Ur = 20 one 
observes more a double diagonal structure with main con¬ 
tributions for positions such that X 2 ^ xi ± 20. 


The energy dependence of both ^e and for all four 
cases of Fig.[T0|and all system sizes between 144 and 1597 
is shown in Fig. 11 As in the Hubbard interaction case 
(see Fig. 1^ the typical values of ^e and increase sys¬ 
tematically with the system size and for each case there 
is a certain narrow, quite well defined, energy band for 
strongly delocalized eigenstates. 


In addition to this, for the three cases presented in the 
three lower rows of panels in Fig.j^one does not see many 
data points for strongly localized states (with ^ 1) in¬ 
side or close to this narrow energy band in contrast to Fig. 

where a lot of eigenstates with very small values of ^e or 
^x are visible (for three out of four cases). The reason for 
this is that the total energy for these three cases is outside 
the interval |A| < 6 for non-interacting product states (at 
A = 2.5) where the two particles are localized more or less 
far away with only small (or absent) effects due to the in¬ 
teraction. Therefore contributions of such products state 
cannot be seen for the particular narrow energy bands 
visible in Fig. El 

In principle this argument also applies to the first row 
of panels in Fig.[^(with Ur = 20 and U = 14.0), i.e. here 
products states with particles localized far away cannot be 
not seen as well. However, for the long interaction range 
Ur = 20 and due to the fact that the interaction is uni¬ 
form in this range there are other products states where 
both particles are localized at a distance smaller than Ur 
which is possible due to the small one-particle localization 
length (. = 4.48 < 20. The spatial structure of these kind 
of product states is not modified by the uniform interac¬ 
tion. Therefore they are strongly localized, but obviously 
the energy eigenvalue of such a short range product range 
is shifted by the mean value of the uniform interaction 
U = 14.0 (with respect to the sum of the two one-particle 
energies) therefore explaining that it is possible to find 
such states for energies close to A « 14. This explains 
also that more complicated effects of the interaction, such 
as the creation of strongly delocalized two-particle states. 



N=1597 

+ 


N=987 

X 


N=610 

!K 

+ 

N=377 

N=233 

□ 

♦ 

,N=144 

O 


U=14.0 

UFt=20,w=0 

bosons 

jjU 







E E 




Fig. 11. Inverse participation ratio of eigenstates versus eigen¬ 
value energy for the system sizes N = 144, 233, 377, 610, 987, 
1597 and the same four cases with Ur > 1 as in Fig. 10 (see 
labels in panels for the values of the parameters U, Ur, w and 
boson or fermion case). The left column of panels correspond 
to the inverse participation ratio in energy representation 
and the right column to the inverse participation ratio G in 
position representation. 


happen if both particles are at an approximate distance 
^ 20 such that the interaction coupling matrix elements 
(between non-interacting product states with both parti¬ 
cles at critical distance ^ Ur) have a more complicated 
and subtle structure due to complicated boundary effects. 
One may note that this particular type of interaction is 
similar to the bag model studied in mM- 

Fig. shows in a double logarithmic scale the size 
dependence of the maximal inverse participation ratios f^E 
(left column) or f^x (right column) for the above and many 
other selected cases, with different values of U, Ur, w 
and boson/fermion case. The typical values of f,E and 
clearly increase strongly with system size N with typical 
exponents b ^ 0.7-1 obtained from the power low fit ^x = 
aN^ as can be seen in Table For two particular cases 
the behavior is even linear with high precision with 6 = 1 
and a fit error below 0.03% (the two data sets shown with 
b = 1.000 ± 0.000 in Table[§. 
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Fig. 12. Largest inverse participation ratio (for a given value 
of N) of FIKS eigenstates versus system size N using all Fi¬ 
bonacci numbers between 55 and 1597 for selected cases of long 
range interactions (same data sets as in Table [^. The left col¬ 
umn corresponds to the inverse participation ratio in energy 
representation and the right column corresponds to the inverse 
participation ratio in position representation. Top (center) 
panels correspond to the boson case with the decay parame¬ 
ter lu = 0 (rc = 1). Bottom panels correspond to the fermion 
case with the decay parameter w = 0. The three numbers in 
the color labels in the top left corner represent the interaction 
strength U, the interaction range Ur and the approximate en¬ 
ergy eigenvalue E (for N = 1597 and the state with largest 
^e)- Note that for a given value of N and set of interaction 
parameters the eigenstates with maximal and maximal ^e 
may be different. 


Actually, these two cases are also characterized by the 
absence of strongly localized states with « 1 in the nar¬ 
row energy band (and accessible by the Arnold! method) 
in a similar way as the case U = 7.2 for Ur = 1 dis¬ 
cussed previously and one may conjecture that the pres¬ 
ence of strongly localized products states (accessible by 
the Arnold! method and with a modest distance between 
both particles) at the same energies as the FIKS eigen¬ 
states might be a necessary condition to lower the expo¬ 
nent from the linear behavior & = 1 to a fractal value 6 < 1, 
eventually due to some weak coupling of FIKS states to 
strongly localized pairs. Such localized pairs with modest 
distance would also be reasonable for the appearance of 
satellite peaks visible in many (but not all) FIKS eigen¬ 
states (see discussion in Sec. [^. 

For certain other cases of Table the exponents are 
clearly below 1, e.g. 6 « 0.7 or 6 « 0.8 indicating a kind 


Table 3. Approximate energy E and results of the power law 
fit ^x = a for selected cases of long range interactions (same 
data sets as in Fig. 121. The bottom six rows of the table 


correspond to the fermion case and the other top rows to the 
boson case. 


u 

Ur 

w 

E 

a 

b 

3.1 

2 

0 

5.72 

0.633 

± 

0.117 

0.942 

± 

0.032 

6.0 

4 

0 

-0.96 

0.637 

± 

0.046 

0.999 

± 

0.013 

7.5 

3 

0 

9.82 

0.465 

± 

0.007 

1.002 

± 

0.003 

8.0 

5 

0 

8.80 

0.753 

± 

0.055 

0.978 

± 

0.013 

10.5 

4 

0 

5.96 

0.684 

± 

0.001 

1.000 

± 

0.000 

14.0 

20 

0 

14.00 

1.118 

± 

0.172 

0.885 

± 

0.027 

8.5 

2 

1 

1.34 

0.727 

± 

0.062 

0.986 

± 

0.015 

17.0 

7 

1 

10.23 

0.339 

± 

0.084 

1.032 

± 

0.043 

18.5 

5 

1 

2.53 

0.485 

± 

0.100 

0.981 

± 

0.036 

22.5 

4 

1 

-1.01 

0.635 

± 

0.001 

1.000 

± 

0.000 

22.5 

5 

1 

4.31 

0.842 

± 

0.168 

0.936 

± 

0.034 

2.0 

6 

0 

1.81 

0.991 

± 

0.147 

0.873 

± 

0.026 

3.5 

2 

0 

5.34 

0.696 

± 

0.111 

0.947 

± 

0.027 

3.6 

6 

0 

2.62 

1.308 

± 

0.158 

0.807 

± 

0.021 

7.0 

6 

0 

2.24 

2.349 

± 

0.718 

0.683 

± 

0.053 

7.5 

5 

0 

7.41 

0.793 

± 

0.138 

0.945 

± 

0.030 

10.9 

5 

0 

11.54 

0.896 

± 

0.141 

0.949 

± 

0.027 


of modest fractal structure of the eigenstates in a similar 
way as for the Hubbard case with U = 7.8 and E « —2.8. 

Furthermore, both the figure labels of Fig, and also 
Table [^provide the approximate energy values for the nar¬ 
row energy delocalization band and in many cases these 
energy values also lie inside the interval \E\ < 6 of non¬ 
interacting product states with both particles localized far 
away, confirming that the strong delocalization effect may 
happen for both cases \E\ < 6 and \E\ > 6. 


7 Momentum and energy representation of 
eigenstates 


It is illustrative to present the FIKS eigenstates which are 
delocalized along the diagonal Xi « X 2 in other represen¬ 
tations such as a momentum representation using discrete 
Fourier transform or in the energy representation in terms 
of non-interacting product one-particle eigenstates, a rep¬ 
resentation already used for the algorithm of the Green 
function Arnold! method described in Sec.l^and Appendix 

m 

We first write a two-particle eigenstate with wave func¬ 
tion 'ip(xi,X 2 ) for a;i, 0:2 S {0, ..., Af — 1} in momentum 
representation by discrete Fourier transform: 

'0(Pi,P2) = ^^ eyi^{ikp^xi+ikp^X2)i’{xi,X2) (13) 

Xi ,X2 


with kp^ = 2'Kpj/N for pj = 0, ..., Af — 1 and j = 1,2. 
The momentum eigenfunction (13) can be efficiently eval¬ 
uated using Fast Fourier Transform using the library fFtw3 
[30] which also works very well with optimal complex¬ 
ity 0{N‘^ log(A^)) (for a two-dimensional discrete Fourier 
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transform) for arbitrary values of N, even for prime num¬ 
bers and not only for powers of two. However, it turns 
out that the density plot of the momentum eigenfunction 
(13) has typically a quite complicated or bizarre structure 
and does not reveal much useful insight in the delocal¬ 
ization effect visible in position representation. Actually, 
the momentum representation with the simple ordering of 
momenta kp with p = 0, ..., N — 1 is not appropriate to 
study the quasiperiodic potential Vi{x) = A cos(aa: -|- (3). 

To understand this more clearly let us revisit the eigen¬ 
value equation of an eigenfunction (j){x) with eigenvalue e 
for the one-particle Hamiltonian with this quasiperiodic 
potential: 

e (f){x) = t[(j){x + 1) + (/'{x — 1)] -I- A cos(Q;a:: + P)(j){x) (14) 

where we have used a generalized hopping matrix element 
t and where for simplicity x may take arbitrary integer 
values for an infinite system and is an irrational 

number such as the golden ratio Q!/(27r) = {pjh — l)/2. 
In Ref. @] a duality transformation was introduced by 
expanding the eigenfunction in the form: 

^{x) = E eyi]>[i{^x apx j3p)\(j){p) (15) 

p 


where the sum runs over all integer values of p, (3 is some 
arbitrary parameter and for convenience we have taken 
out a phase factor exp{if3p) from the precise definition of 
4‘{p). This expansion defines unique coefficients (pip) only 
for irrational values of a. Inserting (15) into (14) one finds 
that the function obeys a similar eigenvalue equation 
of the form: 


e (j){p) = t[(j){p -f 1) -I- 4>{jp - 1)] -I- A cos(q!P -I- P)(t){p) (16) 

with e = 2telX, A = 4t^/A and j3 is the parameter used 
in (15). For |<| = 1 this transformation maps the case 
A > 2 to the case A = 4/A < 2. In Ref. [4], using this 
transformation together with Thouless formula (and some 
technical complications related to a finite size and rational 
approximation limit of a), it was argued that for A > 2 
the eigenfunctions (j){x) are localized with a localization 
length £ = l/log(A/2) and for the dual case (with A < 2) 
the functions (j){p) are delocalized. 

The important lesson we can take from the duality 
transformation (15) is that it uses only a sum over discrete 


momentum values Qp = {/3 + a p) mod (27r), i. e. 




(17) 


relation have to be considered as “neighbor” values in dual 
space, i.e. the natural proper ordering of momentum val¬ 
ues is given by the discrete series qp = {j3 + ap) mod (27r) 
with increasing integer values for p. 

Let us now consider the case of finite system size N 
with periodic boundary conditions ^(0) = (j){N) in ( |l4| ). If 
we want to construct a proper dual transformation for this 
case we have to chose a rational value for a/(2'K) = M/N 
where 0 < M < N and the integer numbers M and N are 
relatively prime (if M and N are not relatively prime we 
would have a periodic potential with a non-trivial period 
being shorter than the system size requiring an analysis 
by Bloch theorem etc.). In this case we may directly use 
(15) to define the duality transformation provided that 
the sum is limited to the finite set p = 0, ..., — 1 [and 

not infinite as for the case of infinite system size with 
irrational Q:/(27r)]. Furthermore, for convenience we chose 
the parameter (3 = 0. Then the discrete momentum values 
Qp become 


Qp = (ap) mod (27r) = 27r 


(pM) mod N 

~n' 


= k. 


v(p) 


(18) 


where kp = 2 ttp/N is the momentu m va lue for the discrete 
Fourier Transform [see also below s and with (t{p) = 
(pM) mod N being a permutation of the set {0, ..., iV—1} 
because M and N are relatively prime. We remind that 
for the eigenstate analysis in the previous Sections we 
had used the choice M = /„_i and N = fn where /„ 
is the n-th Fibonacci number and we note that two sub¬ 
sequent Fibonacci numbers are indeed always relatively 
prime. For this particular choice we call the permutation 
ct(p) the golden permutation. The permutation property 
of crip) Eq. (18) ensure that the discrete momentum 


values Qp of the dual transformation (15) coincide exactly 
with the discrete momentum values used for the discrete 
Fourier Transform for a finite lattice of size N. However, 
there is a modified ordering between qp and kp because 
of the permutation and “neighbor” momenta kp and fcp+i 
of the discrete Fourier Transform are not neighbor val¬ 
ues for the dual transformation and therefore the direct 


naive momentum representation (13) is not appropriate. 


The proper dual transformed representation corresponds 
to the golden permutation Fourier representation defined 

by 


^g(Pl,P2) 


l/(cr(pi),cr(p2)) (19) 

expiiqp^Xi+iqp^X2)ip{xi,X2) 

Xi ,X2 


instead of a continuous integration over q € [0,27r[ which 
would normally be the proper way to perform a Fourier 
transform from the discrete infinite one-dimensional in¬ 
teger lattice space for x to the continuous variable q € 
[0, 27r[. However, the quasiperiodic potential only couples 
(in the dual equation) momenta q and q such that g = (gi 
a) mod (27r) and therefore the discrete sum in (15) is suf¬ 
ficient. Furthermore, two momentum values obeying this 


where t he s econd identity with qp (instead of kp) is valid 
due to (18). For ipgipi,P 2 ) neighbor values in pi or p 2 
correspond indeed to neighbor values in the dual transfor¬ 
mation. 


We mention that for a finite system size N and an irra¬ 
tional choice of a/(27r) the momenta, qp = (ap) mod (27r), 
used for the duality transformation do not coincide exactly 
with the discrete momenta of the discrete Fourier trans- 
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form, in particular the quantity 

' Npa 


t(p) = 


27r 


mod N 


( 20 ) 


would typically not be an integer number. At best one 
could try to define an approximate duality transformation 
with a modified permutation by rounding (201 to the next 


integer number but even in this case one would typically 
not obtain a permutation and it would be necessary to 
correct or modify certain a(p) values in order to avoid 
identical cr(p) values for different integers p. 

If we want to choose a finite system size N which is not 
a Fibonacci number we could try for the choice of a/(27r) 
a rational approximation M/N of the golden ratio (-s/S — 
l)/2 with M being the closest integer to l)/2 and 

the denominator fixed by the given system size. However, 
in this case one might obtain a value of M such that M 
and N are not relatively prime and (if we want to keep the 
same denominator) it would necessary to chose a different 
value of M relatively prime to N and still rather close to 
N{^/5 — 1)/2 therefore reducing the quality of the rational 
approximation. For this reason we have in the preceding 
Sections mostly concentrated on the choice of Fibonacci 
numbers for the system size such that we can use the best 
rational approximation for the golden number and where 
we can always dehne in a simple and clear way the golden 
permutation by a{p) = (p/„_i) mod/„. 

In Fig. the three eigenstates with maximum for 
N = 610, Ur = 1 and the two cases U = 4.5 and U = 7.8 
(and A « —2.8) are shown in the golden permutation 
Fourier representation. One sees clearly that for the cen¬ 
ter of mass coordinate there is a strong momentum local¬ 
ization around a few typical values while for the relative 
coordinate all momentum values seem to contribute to the 
eigenstate leading to momentum delocalization in this di¬ 
rection. This is just dual to the typical behavior of such 
eigenstates in position representation with delocalization 
in the center of mass coordinate and localization in the 
relative coordinate. However, the precise detailed struc¬ 
ture, in momentum space on a length scale of a few pixels 
and well inside the stripes seen in Fig. is still quite 
complicated and subtle. 

The “localization length” in momentum space for the 
center of mass coordinate is considerably shorter for the 
case U = 4.5E with about 10 pixels (i. e. discrete mo¬ 
mentum values) than for the other case U = 7.8 (and 
E « —2.8) with about 30 pixels. This observation relates 
to the stronger quasiperiodic hole-peak structure in the 
eigenstates seen in Figs. |4][^ for the case U = 7.8 (and 
E « -2.8). 

We have also tried for the irrational case and non- 
Fibonacci numbers for N to define an approximate golden 
permutation which in principle provides similar figures as 
in Fig. but with a considerable amount of additional 
irregularities concerning the momentum structure etc. 

Another type of interesting eigenvector representation 
is obtained by an expansion of a two-particle eigenstate in 
the basis of non-interacting one-particle product 

eigenstates. Fig. shows black and white density plots 



Fig. 13. Density plot of the three FIKS eigenstates in golden 
permutation Fourier representation with largest values of 
for N = 610, Ur = !,[/ = 4.5 (left column) or U — 7.8 
(right column). The corresponding energy eigenvalues and val¬ 
ues for both types of inverse participation ratios are: Top 
left: E = -3.09750, Cb = 249.137, = 271.208. Top right: 

E = -2.78586, £,e = 211.058, = 194.241. Center left: 

E = -3.09964, £,e = 239.312, = 265.885. Center right: 

E = -2.78599, ^e = 200.958, = 176.454. Bottom left: 

E = —3.09815, ^e = 233.773, = 250.700. Bottom right: 

E = -2.78593, ^e = 190.171, = 193.885. 


for the amplitudes of certain eigenstates in such a repre¬ 
sentation for the two sizes N = 233 and N = 610 and the 
two values of the interaction U = 4.5 and U = 7.8 (both 
for Ur = 1). Both axis correspond to the one-particle in¬ 
dex ordered with respect to increasing values of the corre¬ 
sponding one-particle energy. We remind that in the sec¬ 
ond variant of the Green function Arnoldi method the 
main calculations are actually done in this energy repre¬ 
sentation, which is therefore more easily accessible than 
the standard position representation. 

One observes a kind of self-similar structure with (ap¬ 
proximate) golden ratio rectangles of different sizes along 
the diagonals. The inverse participation ratio in energy 








16 


K.M.Frahm and D.L.Shepelyansky: Freed by interaction kinetic states in the Harper model 



Fig. 14. Density plot of the FIKS eigenstates in non¬ 
interacting energy representation with the largest value of 
for N = 233 (top panels) or N = 610 (bottom panels), Ur = 1, 
U = 4.5 (left column) or U = 7.8 (right column). Black rep¬ 
resents maximum, grey medium and white minimum values 
of the expansion amplitudes of the shown eigenstate with re¬ 
spect to non-interacting energy product eigenstates 
The horizontal (vertical) axis corresponds to the index i/ (/r) 
ordered with respect to increasing values of the correspond¬ 
ing one-particle energy (c/i) of the first (second) parti¬ 
cle. Top panels for N = 233 correspond to full eigenstates 
and bottom panels for N = 610 correspond to a zoomed re¬ 
gion of size 200 x 200 with left bottom corner at position 
xi = X 2 ~ 100. The corresponding energy eigenvalues and 
values for both types of inverse participation ratios are: Top 
left: E = -3.09669, = 107.409, = 106.818. Top right: 

E = -2.78569, ^e = 117.697, = 102.577. Bottom left: 

E = -3.09750, £,e = 249.137, = 271.208. Bottom right: 

E = -2.78586, £,e = 211.058, = 194.241. 


representation corresponds approximately to the number 
of black dots in the black and white density plots of Fig. 

[H 

We mention that when the one-particle eigenstate or¬ 
dering in the energy representation is done with respect 
to the maximum positions of the one-particle eigenstates 
(instead of the one-particle energy) one obtains a clear 
banded structure with main values/peaks for « /i ± 5 
(Figure not shown). 


8 Implications for cold atom experiments 

Motivated by recent experiments on cold atoms [21] we 
present also some results for a modihed value of the flux 
parameter a used in the quasiperiodic potential Vi{x). In 
the experiment of Ref. [21] the rational value for Q;/(27r) « 




0 5 10 15 20 0 5 10 15 20 

U U 


Fig. 15. Inverse participation ratio ^ipr and variance length 
{S'^x)^^'^ = {{x — of the time evolution two-particle 

state for system size N = 512, iteration time t = 5120 and ai 
(left panel) or 0:2 (right panel) versus interaction strength U. 
The initial state at t = 0 is localized with both particles in 
the center position xo = N/2. Inverse participation ratio and 
variance length have been calculated from an effective one- 
particle density without a center box of size 20% (with respect 
to system size). The data points for 2 < U < 2.5 have been 
calculated with a doubled iteration time t = 10240. The values 
of ai, a 2 from (221, (231 correspond to experimental conditions 
of [^. 


532/738 = 266/369 was used. This value has the finite 
continued fraction expansion [0; 1, 2,1,1, 2,1,1, 8] with 

[00:01,02,03,...] =aoH- - -. ( 21 ) 


02 H-;- 

03 -I- 

We define two numbers aj, j = 1,2 such that ajj{2 tt) is 
irrational and close to the experimental rational value by 

^ = [0; 1,2,1,...] = ^ = 0.7207592200561264... 

ZTT 3 

( 22 ) 

and 

g = [0; 1,2,1,1, 2,1,1,8,...] (23) 

169 

-= 0.7208720926598791... 

where the initial pattern of shown coefficients in the con¬ 
tinued fraction expansion (except the leading zero) repeats 
indefinitely with a period of 3 (or 8) for the case of ai (or 
02 ). The first choice provides a “stronger” irrational num¬ 
ber for Q!i/(27r) while the second choice is closer to the 
experimental value. In this Section we choose for all nu¬ 
merical computations one of these two values (or rational 
approximations of them for the eigenvector calculations) 
and furthermore we hx the phase offset and the interaction 
range by /3 = (-s/S — l)/2 and Ur = 1. 

First, we performed the time-evolution analysis already 
described in Sec. [fusing either ai or 02 . Fig. [Ts] shows 
the dependence of the inverse participation ratio ^ipr and 
the variance length [both computed without the 20% cen¬ 
ter box, see ([^ and @] on the interaction strength U 
{0 < U < 20) Tor a system size N = 512 and an iteration 
time t = 5120. As in Sec. [4] we chose for t = 0 an ini¬ 
tial state with both particles localized at the center point 
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Table 4. Time evolution parameters for the interaction values 
U = 2.25 and U = 3.6 and both values of a using the data sets 
of Fig. |15| These interaction values correspond to the local 
maxima of the squared tail norm ||'i/’taii(t)|p. Note that for U 
close to 2.25 the local maxima visible in Fig. of the length 
scale CipR (computed without the 20% center box) correspond 
actually to 17 = 2.3 with slightly larger values than for U = 
2.25. 


a 

u 

t 

^IPR 

{H) 

S^E 

IIV'tail(t)f 

ai 

2.25 

5120 

32.27 

-4.744 

0.475 

0.00884 

ai 

2.25 

10240 

79.48 

-4.828 

0.127 

0.107 

ai 

3.6 

5120 

101.12 

-0.893 

0.159 

0.0449 

01.2 

2.25 

5120 

30.18 

-4.717 

0.587 

0.00562 

01.2 

2.25 

10240 

64.58 

-4.826 

0.140 

0.0709 

01.2 

3.6 

5120 

45.45 

-0.878 

0.215 

0.0188 


Xq = N/2. For both a values we observe strong peaks for 
both length scales at values U = 2.25-2.3 and U = 3.6 
indicating the possible existence of FIKS states at these 
interaction values (or very close). A closer inspection re¬ 
veals that the first peak close to 17 = 2.25 requires a longer 
iteration time t = 10240 in order to provide saturation of 
the two length scales and therefore in Fig. the data 
points for 2 < 17 < 2.5 are computed with this increased 
iteration time. 

Table 1^ summarizes the results of the quantities ^ipr, 
(H), 6‘^E and ||'0taii(t)|P (see Sec. for the precise def¬ 
inition of them) at the two peak values U = 2.25 and 
U = 3.6. The values of ^ipr in Table for U = 2.25 and 
t = 10240 do actually not exactly correspond to the first 
local maximum visible in Fig. |15| because ^ipr is maximal 
at 17 = 2.3 while the value of 17 = 2.25 corresponds to the 
local maximum of However, detailed eigenvec¬ 

tor calculation for these two interaction values confirm 
that globally the value U = 2.25 is slightly more optimal 
than U = 2.3 with stronger delocalization. 

In Table 1^ we provide for the case U = 2.25 also the 
results for the two iteration times t = 5120 and t = 10240. 
Obviously, ^ipr and ||V’taii(t)|P are considerably increased 
at t = 10240 but already at t = 5120 the strong delocaliza¬ 
tion FIKS effect is visible. The average energy value of the 
tail state is rather sharp with a modest variance S^E for 
all cases, but also with an additional significant decrease 
of S'^E between t = 5120 and t = 10240 (for U = 2.25). 

Globally Fig.and Tablej^show that the FIKS effect 
is stronger for U = 2.25 but at this value it requires a 
longer iteration time to be clearly visible. This observation 
is also confirmed by Fig. which shows for ai and both 
interaction values U = 2.25 and U = 3.6 the density plots 
and the one-particle density of three time evolution states 
at t = 100, t = 1000 and t = 10000. In both cases the state 
is clearly localized at the beginning at f = 100 and it is 
delocalized over the full system size at t = 10000 (with a 
small weight and along the diagonal Xi Ki x 2 as discussed 
in Sec. |^. However, for the intermediate time t = 1000 
the state for U = 2.25 is considerably less delocalized than 
the state for U = 3.6 at the same iteration time clearly 



Fig. 16. Density plot (three top rows of panels) of time evolu¬ 
tion two-particle states for system size N = 512, the case oi, 
interaction range Ur = 1, interaction strength U = 2.25 (left 
column) or U = 3.6 (right column), iteration times t = 100 
(first row), t = 1000 (second row) and t — 10000 (third row); 
panels show the whole system range (0 < x^X 2 < 512). The 
fourth row of panels shows the one-particle density pi(a;) in a 
semi-logarithmic representation for the same states as in the 
three top rows of panels. The initial state at t = 0 is localized 
with both particles in the center position xq = N/2. 
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Fig. 17. Density plot for the time dependence of one-particle 
density from the time evolution state with a:-position (0 < a; < 
512) corresponding to the horizontal axis and time t {0 < t < 
10240) corresponding to the vertical axis. Here U = 2.25 (left 
panel), U = 3.6 (right panel) and q = ai, Ur — 1. 


confirming the slower delocalization speed for U = 2.25. 
Thus the velocity of FIKS pairs is smaller at 17 = 2.25 than 
at U = 3.6 but the weight of FIKS pairs in the initial state 
is larger at 17 = 2.25. Apart from this the delocalized tails 
of the state at t = 10000 appear somewhat “stronger” 
or “thicker” for U = 2.25 explaining the larger values 
of CiPR (for ai). Note that Fig. [I^ shows the full time 
evolution states while Fig. in Sec. with the golden 
ratio value for a/(27r), shows only a zoomed range for the 
right delocalized branche between the right border of the 
20% center box and the right border of the full system. 

Fig. [IT] shows the time evolution of the one-particle 
density (for ai) with the x-dependence corresponding to 
the horizontal axis and with the t-dependence (0 < t < 
10240) corresponding to the vertical axis. This figure pro¬ 
vides clear and additional confirmation that the delocal¬ 
ization effect is stronger and slower for U = 2.25 than 
for U = 3.6. It also confirms the linear (ballistic) increase 
of the delocalized part of the state with time (see also 
Fig.^. We mention that the other value a 2 provides very 
similar figures as Figs. [T5| and [TB] with a slightly reduced 
delocalization effect for both interaction values. 

Following the procedure described in the beginning of 
Sec.j^we have also computed eigenstates using the Arnoldi 
Green function method with the average energy values 
(H) of Table as initial Green’s function energy for the 
smallest system size. The Green function energies are re¬ 
fined for larger system sizes using the energy eigenvalue 
of a well delocalized eigenstate of the last smaller sys¬ 
tem size. Following the spir it of t he p revious explications 
[see text between Eqs. (17) and (18)] we choose rational 
approximations of ai/ (27r) and a^TI^Tr) using their con¬ 
tinued fraction expansions (22) and (23) which provide 
suitable system sizes given as the denomators of the ra¬ 
tional approximations. Using a minimal (maximal) system 
system size ^ 40 (~ 10000) this provides for ai the values 
N = 43, 111, 154, 265, 684, 949, 1633, 4215, 5848, 10063 
and for a 2 the values N = 43, 369, 412, 1193, 1605, 2798, 
7201, 9999. Note that the system size 369 corresponds 
to the rational approximation a2/(2Tr) « 266/369 used 
in the experiments of Ref. [2T]. For each system size we 



Fig. 18. Density plot of FIKS eigenstates for rational approx¬ 
imations of Q2/(27r) and U = 2.25 (left column), U = 3.6 
(right column), N = 369 (top panels), N = 1605 (center and 
bottom panels); Ur = 1. The corresponding energy eigenval¬ 
ues and values for both types of inverse participation ratios 
are: Top left: E = -4.85051, Ce = 98.462, = 118.308. Top 

right: E = -0.92196, = 113.232, G = 108.389. Center left: 

E = -4.84994, = 428.375, G = 566.237. Center right: 

E = —0.92198, = 309.040, G = 260.125. Bottom panels 

show a zoomed region of size 100 x 100 with left bottom corner 
at position Xi = X 2 = 350 of the center panels. 


use the corresponding rational approximation of aj/(27r) 
(j = 1, 2) and /3 = {'/b — l)/2) to determine numerically 
certain eigenstates by the Green function Arnoldi method. 

In Fig. ^ we show selected strongly delocalized eigen¬ 
states for a 2 and the two interaction values U = 2.25 and 
U = 3.6 and the system sizes N = 369 and N = 1605. All 
eigenstates provide nice FIKS pairs with a quite specific 
particular pattern on the diagonal xi « X 2 which corre- 
ponds, for each of the two interaction values, rather well to 
the pattern of (the delocalized tails of) the time evolution 
states for t = 10000 visible in Fig. For N = 1605 the 
pattern for U = 2.25 seems be to considerably more com¬ 
pact than the pattern for U = 3.6 which is also confirmed 
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Fig. 19. Largest inverse participation ratio (for given values of 
N and approximate energy) of FIKS eigenstates versus system 
size A'’ in a double logarithmic scale for Ur = 1, for rational ap¬ 
proximations of oi / {2tt) (top panels) or a ?2 / {2n) (bottom pan¬ 
els) and for U = 2.25, E ~ —4.85 (left column) or U = 3.36, 
E ~ —0.92 (right column). The used system sizes are 43, 111, 
154, 265, 684, 949, 1633, 4215, 5848, 10063 given by the de¬ 
nominators of the rational approximations of ai/(27r) and 43, 
369, 412, 1193, 1605, 2798, 7201, 9999 for the rational approx¬ 
imations for a2/(27r). The blue line with stars corresponds to 
the inverse participation ratio in position representation, 
the red line with crosses to the inverse participation ratio 
in energy representation and the green line to the power law fit 
^x = aN’’ with fit values given in Table Note that for given 
values of N and approximate energy the eigenstates with max¬ 
imal ^x and maximal may be different. 


Table 5. Results of the power law fit = a for the four 
cases of Fig. |19| 


u 

Ur 

a 

E 

a 

b 

2.25 

1 

ai 

-4.85 

0.424 

± 

0.013 

0.990 

± 

0.004 

3.6 

1 

ai 

-0.92 

1.878 

± 

0.348 

0.685 

± 

0.027 

2.25 

1 

(y.2 

-4.85 

0.418 

± 

0.055 

0.971 

± 

0.018 

3.6 

1 

(y.2 

-0.92 

3.333 

± 

1.197 

0.581 

± 

0.050 


by a considerably larger value of ^x- The eigenstates for 
the ai case are very similar for comparable system sizes. 

Fig. [T^ shows the size dependence of ^x and for the 
four cases corresponding to any combination of the two 
interaction and the two flux values. The fit results of the 
power law fit = a are shown in Table[^ For U = 2.25 
both fits for the two flux values are very accurate with 
exponents 6 « 1. For U = 3.6 the fit quality is somewhat 
reduced and the exponents are quite smaller b « 0.7 for 
ai and b « 0.6 for a 2 indicating a certain fractal structure 
of eigenstates. At N 10000 the maximal values of for 
U = 2.25 at both flux values are at least four times larger 
than the maximal values of ^x for U = 3.6. We also observe 
that the density of good FIKS pairs for [/ = 2.25 and both 
flux values is extremely high. In Secs. for the rational 
approximation of the golden ratio for a/(27r) only the case 



Fig. 20. Density plot of two selected eigenstates for U = 4.5, 
Ur = 1 and for the rational approximations of ai j(2ir') . The 
corresponding system sizes, energy eigenvalues and values for 
both types of inverse participation ratios are: Lejt: N = 369, 
E = -2.21758, ^E = 11.378, ^x = 21.029 (2nd largest value of 
^E and largest value of ^x for this system size and approximate 
energy). Right: N = 1605, E = -2.21949, ^e = 16.367, ^x = 
26.748 (largest value of ^e for this system size and approximate 
energy). The left panel shows the full state of size 369 x 369 
and the right panel a zoomed region of size 369 x 369 with left 
bottom corner at position a:i = 3:2 = 0 and outside the zoomed 
range no data points different from blue (for zero amplitude) 
are visible. 


for U = 7.2 has a comparable density of good FIKS pairs 
(see bottom panels of Fig. [^. 

We have also tested (for 02 ) the interaction strength 
1/ = 4.5 with approximate energy E = —3.1 which pro¬ 
vided nice FIKS pairs for the golden ratio case studied 
in Sec. However, here we should not expect delocal¬ 
ized FIKS pairs since according to Fig. [I^ the value of 
CiPR obtained from the time evolution state is very small. 
On the other side, the variance length shows some mod¬ 
estly increased values and it might be useful to verify such 
cases as well. We applied the standard procedure of en¬ 
ergy refinement with the Green function Arnoldi method 
on t/ = 4.5 with the initial energy E = —3.1 which imme- 
diatedly selected E « —2.2 as “optimal” energy range (to 
maximize ^e)- Despite some modestly delocalized eigen¬ 
states with ^E ~ 15 and ^x 25 (for the largest consid¬ 
ered systems sizes N = 412, 1193 and 1605) there are no 
FIKS pairs with strong delocalization along the diagonal. 
Fig. [20| shows for 02 and N = 369 or TV = 1605 two such 
modestly delocalized eigenstates which have some “cigar” 
form but with a rather short length ^ 50-80 and a rather 
elevated width ^ 20-30. It seems that the variance length, 
in contrast to ■Cipr? does not really allow to distinguish 
between these kind of states and nice FIKS eigenstates. 
Furthermore this example shows that suitable parameters 
U and E for FIKS states depend strongly on the flux pa¬ 
rameter a, an issue which is more systematically studied 
in the next Section. 


9 Dependence on flux values 

A problem with a systematic study of the dependence of 
the FIKS effect on different flux values is to select a suit- 
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able set of irrational numbers of comparable quality and 
which have roughly the same distance. For this we con¬ 
sider at hrst rational numbers p/89 with 44 < p < 88 
where the denominator 89 has the nice feature of being 
both a prime and a Fibonacci number. We compute for 
each of these rational numbers the canonical variant of 
its finite continued fraction expansion |32j , reduce the last 
coefficient by 1 and add an inhnite sequence of entries of 
1. This provides the inhnite continued fraction expansion 
of an irrational number which is rather close to the initial 
fraction p/89 and which has “a golden tail” for the con¬ 
tinued fraction expansion. It turns out that for each value 
of p the difference between p/89 and the corresponding 
irrational number is approximately 5 x 10“® therefore pro¬ 
viding a nice data set of irrational numbers between 0.5 
and 1. 


In particular for p = 55, where 55/89 is a rational 
aproximation of the golden number, we have 55/89 = 
[0; 1,1,1,1,1,1,1,1, 2]. The procedure reduces the last co¬ 
efficient from 2 to 1 and adds the inhnite sequence of unit 
entries just providing exactly the continued fraction ex¬ 
pansion of the golden number (with all coefficients being 
unity). The golden number is therefore one of the data 
points in the selected set of irrational numbers. For p = 64 
we hnd the irrational value 0.7191011235955056 ... which 
is by construction very close to 64/89 but also rather close 
to 266/369 « 0.72087, which was used in the experiment 
of R ef. 121| . and also to the two irrational numbers (221 
and (23) used in the previous Section. 

Using these irrationals values for a/(27r) and (3 = (-s/S— 
l)/2 we have performed the time evolution analysis de¬ 
scribed in Sec. |4 for system size N = 512, iteration time 
t = 5120 and the interaction interval 0 < U < 10 in 
steps of AU — 0.25 providing in total 45 x 21 data sets. 
The main results of this analysis are shown in Fig. 
containing two density plots in a/{2Tr)-U plane for the 
squared tailed norm ||t/’taii(i)|P and the inverse participa¬ 
tion ratio ^ipR (without 20% center box) both providing 
the most reliable measure of delocalization in the frame¬ 
work of the time evolution analysis (the variance length 
provides a considerable amount of fluctuation peaks also 
when the other two quantities are very small as can be 
seen in Figs. and 15). 

Concerning the density plots of a quantity p we men¬ 
tion that we apply the attribution of the different color 
codes to uniform slices of p” with r < 1 being some expo¬ 
nent, of typical choice 1/4 or sometimes 1/8, to increase 
the visibility of small values of p. In Fig. we used for the 
density plot of the squared tail norm the standard choice 
r = 1/4 due to the large ratio ~ 10^^ between maximum 
and minimum values but for ^ipr where this ratio is ~ 10^ 
we chose exceptionnally r = 1. For these plot parameters 
the density plots for these two quantities provide rather 
coherent and similar results for parameter regions with 
strong delocalization. All raw data of Fig. 21 are available 
for download at 1^. 


The density plots of Fig. show that for values of 
a/ (27r) close to the simple fractions 1/2, 2/3, 3/4 and even 
4/5 there is a certain rather uniform delocalization effect 



0.5 0.6 0.7 0.8 0.9 
ct/(2jt) 

8 -_K- 60 


0.5 0.6 0.7 0.8 0.' 

a/{2jt) 



0.6 0.7 

a/(2jt) 


0.7 0.8 0.9 

a/(23i) 


Fig. 21. Top panels : Density plot of the squared tail norm 
||V’taii(t)|P (left panel) or the inverse participation ratio com¬ 
puted without 20% center box (right panel) with horizontal 
axis representing the parameter a/(27r) and vertical axis rep¬ 
resenting the interaction strength U using a time evolution 
state for system size N = 512, iteration time t = 5120, inter¬ 
action range Ur = 1 and a localized initial state for t = 0 with 
both particles in the center position xo = N/2. The bottom 
left panel shows a zoomed range with 0.6 < a/{2n) < 0.75 
and 2 < U < 8 of the top left panel. The two arrows indicate 
the value of the golden ratio Q/(27r) = (\/5 — l)/2 « 0.618 
and the value a/(27r) = 266/369 « 0.721 used in the experi¬ 
ments of Ref. m- The bottom right panel shows the inverse 
participation ratio for U = 0 versus the parameter a/{2n) and 
computed with (red crosses) and without (blue squares) the 
20% center box. The maximal value is for the squared tail 
norm 0.12519 (for a ~ 0.596 and U = 9.5) and for the the 
inverse participation ratio without 20% center box 188.68 (for 
a « 0.753 and U = 4.5). The data of this figure are obtained 
with /3 = (\/5 — l)/2 (and not /3 = 0 as the data of Fig. [^and 
Table in Sec.[^. 


for nearly all interaction values U > 0. We attribute this 
observation to a strong enhancement of the one particle 
location length even in absence of interaction for these flux 
values as can be seen in the bottom right panel of Fig. 
which compares the two variants of the inverse partici¬ 
pation ratio computed with or without the 20% center 
box for vanishing interaction strength U = 0. The first 
variant of ^ipr measures rather directly the effective one- 
particle localization length and is quite enhanced for the 
above simple fractions if compared to the standard value 
£ = l/log(A) « 4.48 for A = 2.5 (for irrational values of 
a 1(211) and infinite system size) [3]. It seems that for the 
irrational values close to simple fractions the system size 
iV = 512 is still too small to see this standard value and 
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one observes an effective enhanced one-particle localiza¬ 
tion length. We have verified this also by direct diagonal- 
ization for some example cases. 

Apart from the simple fractions there are certain com¬ 
binations of a/(27r) and U with a strong FIKS effect and 
a non-enhanced one-particle localization length. For ex- 
emple for the golden ratio case one recovers the peaks at 
U = 4.5 and U = 7.25 (being close to 7.2 found in Sec. 
Q and also for a/(27r) close to the value of 266/369 of 


ef. [H] there are two modest peaks of green color at 
U = 2.25 and U = 3.5 (being close to 3.6 found in the 
previous Section) as can be seen from the zoomed den¬ 
sity plot of the squared tail norm (bottom left panel in 
Fig. 15). We remind that the value U = 2.25 also required 
longer iteration times {t = 10240 instead of < = 5120) to 
be more clearly visible thus explaining the green (instead 
of red) color for this data point since in Fig. 15 we have 
t = 5120. 

Other examples are a « 0.596 and U = 9.5 (with 
maximal value of the squared tail norm of all data sets), 
a ~ 0.753 and U = 4.5 (with maximal value of ^ipr 
without 20% center box) and ol « 0.697 with two in¬ 
teraction values U = 1.75 and U = 2.5. We also com¬ 
puted some eigenvectors by the Green function Arnoldi 
method for these four cases which clearly confirms the ex¬ 
istence of FIKS eigenstates in each case. For example the 
strongest delocalized eigenstate for a « 0.596, U = 9.5 
and N = 1533 corresponds to E = 4.72729, = 426.076, 

= 324.511 and for a « 0.753, U = 4.5 and N = 1837 
to F; = -0.68824, = 3618.270, = 955.650. 

We mention that, for the golden ratio value, the data 
of Fig. 21 are not perfectly identical/coherent to the data 
of Fig. 1^ and Tablej^due to the different phase offset jS = Q 
used for the latter. 


10 Discussion 


The results presented in this work clearly show the ap¬ 
pearance of completely delocalized FIKS pairs induced by 
interaction in the non-interacting localized phase of the 
Harper model when all one-particle eigenstates are expo¬ 
nentially localized. The number of sites (states) ^ popu¬ 
lated by FIKS pairs grows with the system size approxi¬ 
mately like a power law ^ oc with the exponent being 
approximately in the range 0.7 < b < 1. We assume that 
the actual value of b may depend on the energy range and 
interaction strength. It is possible that for 5 < 1 we have 
some multi-fractal structure of FIKS eigenstates. In spite 
of a significant numerical progress and large system sizes 
studied here (we note that the total Hilbert space of the 
TIP problem is Nh = ~ 10® at maximal N = 10946) 

there are still many open aspects in this interesting prob¬ 
lem of interplay of interactions, localization and quasiperi¬ 
odicity. Below we list the main of them. 

Physical origin of FIKS pairs. We see rather subtle 
and complex conditions for appearance of FIKS pairs. 
Their regions of existence are rather narrow on the en¬ 
ergy interval, flux and in the range of interactions (see 
e.g. Figs. |8|21 ). However, at optimal parameters we may 


have up to 12% of states from the initial configuration 
with particles on the same or nearby site being projected 
on FIKS pairs. Thus the optimal conditions and the phys¬ 
ical understanding of the FIKS effect should be clarified. 
If the energy eigenvalue equation of the original Hamil¬ 
tonian Q-® is rewritten in the basis of non-interacting 
eigenstates then it gets the form [5] 


(Cmi + em2)Ami,m2 + ^ Qmi .ma ,m( .m( .m! 

,7712 


(24) 


where Xmi.ma are eigenfunctions of the TIP problem in 
the basis of the non-interacting product states |(/mi, 4 >m 2 > 
introduced in Appendix [B| Note that the (second variant) 
of the Green function Arnoldi method computes rather 
directly Xmi.ma and that f^E is the inverse participation 
ratio in this energy representation. The transition matrix 
elements produced by the interaction are (for the Hubbard 
interaction case) 


Q 


mi,m 2 ,771i ,7712 


XI *^^2 (^) (^) ( 2 ^) ( 25 ) 


with (prnix) =<x\4>rn > being the one-particle eigenfunc¬ 
tions of (IH with the one-particle energies Cm- 

We know that one-particle energies of the Harper model 
at A > 2 have gaps and localized eigenstates. We can as¬ 
sume that the sum of TIP energies also has gaps (or quasi¬ 
gaps) and thus there are some narrow FIKS bands with 
TIP energy width Aeg. On the other side the interaction 
generates some transition matrix elements between these 
band states with a certain typical transition amplitude 
teff oc U. Since the energy inside the FIKS band oscil¬ 
lates quasiperiodically with the distance along the lattice 
we can have approximately the situation of the original 
Aubry-Andre model so that the delocalization transition 
will take place as soon as Agg < 2teg- We think that this 
is the physical mechanism of TIP delocalization in the 
Harper model. However, the concrete verification of this 
mechanism is not so simple: the matrix elements are also 
oscillating with the lattice distance and there are quite a 
several of them (and not only two as in the Harper model), 
there are also energy shifts produced by interaction (the 
diagonal terms) and probably these shifts are at the ori¬ 
gin of narrow regions of interaction where the FIKS pairs 
appear. 

There are some indications from the kicked Harper 
model |33l|34l|35j[36] , that coupling transitions between a 
large number of sites leads to new effects and even ballis¬ 
tic delocalized states. Such ballistic states appear in the 
regime when the classical dynamics is chaotic and diffusive 
and from the analogy with the quantum Ghirikov standard 
map m one would expect to find only pure point spec¬ 
trum of exponentially localized sates. Indeed, there are 
only two transition elements between sites in the Harper 
model while in the kicked Harper model there are sev¬ 
eral of them. The results presented here also indicate that 
the interactions with a longer range have a larger fraction 
of FIKS pairs. Thus for Ur = 5, which has an optimal 
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interaction range, comparable with the one-particle local¬ 
ization length, we obtain a rather large weight of FIKS 
pairs of about 10% in energy and 1 0% in the interaction 
range 0 < ?7 < 20 (see Figs. ggig. These fractions ex¬ 
ceed significantly the typical interaction and energy ranges 
for FIKS pairs with the Hubbard interaction. 


We assume that the spectrum of FIKS pairs has a 
structure similar to the spectrum of the delocalized phase 
in the Aubry-Andre model at A < 2, being close to the 
ballistic spectrum. Indeed, in the time evolution of wave 
packet (see e.g. Fig. we see the lines with a constant 
slope corresponding to a ballistic propagation with a con¬ 
stant velocity. The maximal velocity is Vp « ccmax/imin ~ 
0.2 being smaller then the maximal velocity Vp = 1 for one 
particle at A = 0. It is clear that much more further work 
should be done to obtain a deeper physical understanding 
of the FIKS effect in the Harper model. 


Mathematical aspects. The question about the exact 
spectral structure of FIKS pairs is difficult to answer only 
on the basis of numerical simulations since the system 
size remains always finite and subtle fractal properties of 
the spectrum require more rigorous treatment. There are 
significant mathematical advancements in the analysis of 
quasiperiodic Schrodinger operators reported in [S1I3H1EH] • 
We hope that the results presented here will stimulate 
mathematicians to the analysis of properties of the FIKS 
phase. 


FIKS pairs in cold atom experiments. The results pre¬ 
sented in Sec. show that the FIKS pairs exists at the 
irrational flux value a/(27r) « 532/738 realized in the re¬ 
cent experiments |21j . However, the initial state prepared 
in |2Ij had approximately one atom per each second site 
thus being rather far from the initial configuration con¬ 
sidered here. We think that an initial state with all atoms 
located in the center of the lattice will be much more fa¬ 
vorable for the observation of FIKS pairs. Indeed, such a 
state is rather similar to the initial state considered in our 
paper (two particles on same on nearby sites) and thus we 
expect that in the experiments one will see ballistic prop¬ 
agating FIKS pairs on the tails of probability distribution 
like it is well seen in Figs. 

We note that the initial state with all atoms in the cen¬ 
ter of the lattice had been used in cold atoms experiments 
in the regime of the Aubry-Andre model [50] . In these ex¬ 
periments a subdiffusive delocalization of wave packet has 
been observed being similar to the numerical studies of the 
nonlinear Schrodinger equation on the disordered lattice. 
Indeed, in the center of the packet with many atoms the 
Gross-Pitaevskii description can be more adequate com¬ 
paring to the TIP case considered here. However, on the 
tails of probability distribution on larger distances from 
the center there are only a few atoms and only FIKS pairs 
can reach such far away distances. Thus it is rather pos¬ 
sible that the probability tails will contain mainly FIKS 
pairs. In fact the experimental data in [20] (Fig. 3a there) 
have a plateau of probability at large distances. However, 
at present it is not clear if this is an effect of fluctuations 
and experimental imperfections or a hidden effect of FIKS 
pairs. We think that the present techniques of experiments 


with cold atoms in quasiperiodic lattices allow to detect 
experimentally the FIKS pairs discussed in this work. 

FIKS pairs for charge-density wave and high mate¬ 
rials. We can expect that at finite electron density in a ID 
potential at certain conditions the main part of electrons 
below the Fermi energy will remain well localized creat¬ 
ing an incommensurate quasiperiodic potential for a small 
fraction of electrons in a vicinity of the Fermi level. The 
FIKS pairs can emerge for this fraction of electrons. Such 
situations can appear in the regime of charge-density wave 
in organic superconductors and conductors at incommen¬ 
surate electron density created by doping (see e.g. [40]). In 
such a regime it is possible that the FIKS pairs will give a 
significant contribution to conductivity in such materials. 
The proximity between the charge-density wave regime 
and high T(. superconductivity in cuprates BDiia also in¬ 
dicates a possibility that FIKS pairs can play a role in 
these systems. However, a more detailed analysis of finite 
density systems is required for the solid state systems. 

We think that the various aspects of possible implica¬ 
tions of FIKS pairs in various mathematical and physicals 
problems demonstrate the importance of further investi¬ 
gations of this striking phenomenon. 

This work was granted access to the HPC resources of 
CALMIP (Toulouse) under the allocation 2015-P0110. 


A Description of the Arnoldi method 

For both Lanezos and Arnoldi methods one chooses some 
initial vector |/'i >, which should ideally contain many 
eigenvector contributions, and determines a set of orthonor¬ 
mal vectors |Ci > • ■ •; ICn .4 >j where we call n^. the Arnoldi 
dimension, using Gram-Schmidt orthogonalization on the 
vector II\fk> with respect to |Ci > ..., |Cfc > to obtain 
ICfc+i >• This scheme has to be done for fc = 1, ..., 
and it also provides an approximate representation ma¬ 
trix of “modest” size ha x ua oi H on the Krylov sub¬ 
space generated by these vectors. The largest eigenvalues 
of this representation matrix, also called Ritz eigenvalues, 
are typically very accurate approximate approximations of 
the largest eigenvalues of H and the method also allows 
to determine (approximate) eigenvectors. It requires that 
the product of H to an arbitrary vector can be computed 
efficiently, typically for sparse matrices H but, as we will 
see in the next Section, even non-sparse matrices such as 
resolvent operators can be used provided an efficient algo¬ 
rithm for the matrix vector product is available. 

In its basic variant the Arnoldi method provides only 
the eigenvalues and eigenvectors for the largest energies 
(in module) at the boundary of the band which is not at 
all interesting and in our case it is indeed necessary to 
be able to determine accurately the eigenvalues close to a 
given arbitrary energy. 

The standard method to determine numerically a mod¬ 
est number of eigenvalues localized in a certain arbitrary 
but small region of the eigenvalue space for generic large 
sparse matrices is the implicitly restarted Arnoldi method. 
In this method the initial vector is iteratively refined by 
removing eigenvector contributions whose eigenvalues are 
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outside the energy interval of interest using a subtle proce¬ 
dure based on shifted QR-steps [53]. Using this algorithm 
we have been able to determine eigenvalues and eigenvec¬ 
tors for system sizes up to = 700-1000 but the compu¬ 
tation time is very considerable due to the large number 
of iterations to achieve convergence of eigenvectors. Fur¬ 
thermore, in order to limit the computational time to a 
reasonable amount one has to accept eigenvalues of mod¬ 
est quality with = 10“^^-10“® where the quantity 

5^E{il;)=<'il)\[E-Hf\'il)> (26) 

measures the quality of an approximate eigenvector |'(/i > 
with an approximate eigenvalue E =<ip\ iJ |'0>- Writing 

>= -\-s\6'tp ^ with 771"0exact 7/exact | V^exact ^ 

and ||(5^/’|| = 1 one finds easily that < V'exactl^V^ >= 

(due to normalization of |'0> and |'i/'exact>) and therefore 
E = inexact + O(e^) and S'^E^ip) = 0{e^). Therefore a 
value of 5‘^E[ip) = 10“® implies e ~ 10“"^. 


B Details of the Green function Arnoldi 
method 


In this appendix we provide some of the details concern¬ 
ing the Green function Arnoldi method. For the sake of 
simplicity, we will omit (most of) the details concerning 
the (anti-)symmetrization of two-particle states for bosons 
(fermions) and the corresponding matrix operators act¬ 
ing on them. These details are of course important and 
must be dealt with care and precision when implement¬ 
ing the algorithm. For example, the efficient algorithm for 
the position-energy transformation (see below) requires a 
temporary extension of (anti-)symmetrized states of the 
boson (fermion) space of dimension N 2 to states in the 
general non-symmetrized two-particle space of dimension 
and a corresponding reduction afterwards. However, 
the details for this kind of extensions or reductions with 
eventual y/2 factors etc. are based on the application of 
basic text book quantum mechanics and would only ob¬ 
scure the following description. 

Our algorithm exploits the fact that the interaction 
operator U acts only on a small number of sites {2U]i — 
1) N N'^ |57| given by the set 


S = 



\xi 



(see again [55]). Let us denote by 


P= X! \xi,X2><Xi,X2\. 

{Xi ,X2)GS 


(27) 


(28) 


the projector on the sites belonging to the set S. Obviou^ 
P commutes with the interaction operator U given in M 
and we have PUP = PU = UP = U. For the case of 
the Hubbard interaction with Ur = 1 we even have U = 
UP where U is the interaction strength and corresponds 
to the situation considered in [TH1IT7| . However for Ur > 
1 and re > 0 we note that the operators U and P are 


not proportional (but of course they still commute). We 
denote by Hq = the Hamiltonian in absence 

interaction and by Gq = {E — Hq)~^ the Green function 
or resolvent of Hq. Furthermore we denote by Gq = PGqP 
the projected resolvent (for U = 0) which is a non-trivial 
(non-zero) operator only with respect to its diagonal block 
associated to the subspace corresponding to the set S. 

In this case we can state the following “magic” exact 
formula (§ which is the basic ingredient of our numerical 
approach. This formula can be obtained from a perturba¬ 
tive expansion of G with the interaction as perturbation 
and an exact resummation of all terms except the first 
one. It is also possible to provide an algebraic direct proof 
without use of an expansion and we insist on the fact that 
Q is exact and not approximate. Details for both deriva¬ 
tions are given in Appendix [C] 

The key for an efficient determination of G\(p> using 
(j^ is the observation that the operator (1 — UGo)~^U 
applied to any vector provides only non-zero contribu¬ 
tions on the subspace associated to the set S and the 
matrix inverse is done for a matrix of size UrN <C N 2 
[or (Ur — 1)N N 2 for the fermion case] [5^ once Go 
has been determined. This approach generalizes an idea 
already used in [mim where (for the case of Hubbard 
interaction) the projected resolvent (for arbitrary U) G = 
PGP = Go (1 — t/Go)“^ was calculated to determine the 
localization properties of two interacting particles in one 
dimension from G (we remind that in |16Lll7j a disorder 
and not quasiperiodic potential was studied). 

The numerical algorithm to determine efficiently G\ip> 
is composed of two parts. The first part is to calculate Go 
and the matrix inverse (1 — UGo)~^ which needs to be 
done only once if the value of E is not changed. The sec¬ 
ond part is to evaluate efficiently the successive matrix 
vector products (with Gq, G, (1 — GGo)“^ etc.) accord¬ 
ingly to the formula 

For both parts we need first to diagonalize the one- 
particle Hamiltonian h resulting in eigenvectors > and 
eigenvalues ei, which can be done with complexity 0{N^) 
(or even better using inverse vector iteration for the eigen¬ 
vectors). Then the resolvant Gq can be determined from 


<a;i,X2| Go 


|2/i,J/2>= 


((lyjxi) <(>^( 3 ^ 2 ) (PtJ.{y2) ((ujyi) 
E - e„ - ef_t 


g{E]x,y) 


Y - e ^; X2, 92) 


E 


E - 


=<x\{E-h) ^\y> 


(29) 

(30) 


where g{E; x, y) is the one-particle Green function and 
4>v{x) =<x\l^u>- 

We use (29) to determine the projected resolvent Go, 
i. e. for (xi, X 2 ), ( 2 / 1 , 2 / 2 ) G S. This requires only 0{N^ G|.) 
operations in total since for each value of v we can deter¬ 
mine the one-particle Green function as inverse of a tridi¬ 
agonal matrix (with periodic boundary conditions) with 
0{N‘^) operations using a smart formulation of Gauss al¬ 
gorithm. Then, still for the same value of v, we have to up- 








24 


K.M.Frahm and D.L.Shepelyansky: Freed by interaction kinetic states in the Harper model 


date the sums for all possible values {xi,X 2 ), (yi, 2 / 2 ) G S 
which costs 0{N‘^ C/|.) operations which is dominant (or 
comparable if Ur = 1) to the complexity of the one- 
particle Green function evaluation. The sum/loop over 
V leads then to a further factor of N giving 0{N^ C/^) 
operations. The subsequent matrix inverse to determine 
(1 — UGo)~^ requires 0{N^ U^) operations. We mention 
that for the Hubbard interaction case Ur = 1 this algo¬ 
rithm to determine Gq and the inverse was already imple¬ 
mented and explained in Ref. HZ]. 

For the second part of the algorithm we still need an 
efficient method to evaluate Go|(/5 > for a given vector 
\lp>. This can actually be done by a transformation from 
position to energy representation, i. e. an expansion of 
\ip> using the eigenvectors of Hq given as product states 
\4’v-,4’tJ. >• This transformation can be done with essen¬ 
tially 0{N^) operations using the trick to transform first 
the coordinate of the first particle and then in a sepa¬ 
rate subsequent step the coordinate of the second particle. 
Each one-particle transformation requires 0{N‘^) opera¬ 
tions but it has to be done for N possible positions of the 
other particle and the transformation for the other particle 
gives a further factor of 2 resulting in ~ 2N^ addition and 
multiplication operations for one two-particle transforma¬ 
tion. The transformation back into position representation 
can be done similarly. 

Since Go is diagonal in the energy representation (with 
eigenvalues (E — — e^)“^) the product Go|(/3 > in this 

representation only requires G(iV^) operations. Once this 
is done the resulting vector is transformed back into po¬ 
sition representation (also with 0{N^) operations). Then 
the product of the matrix (1 — UGq)~^U to a vector in 
position representation only requires 0{N‘^U'^) operations 
(provided that the matrix inverse is calculated and stored 
only once in advance for a fixed value of E). Finally a 
further double-transformation-multiplication step with Gq 
is necessarv Combing all this it is possible to evaluate 
G\ip> by ^ by 0{N^) operations (but with a rather big 
prefactor) where the most complex part consists of the 
two position-energy transformations and the two inverse 
energy-position transformations. 

In summary we have described an algorithm to deter¬ 
mine G\(p> by 0{N^ U^) operations for the initial prepa¬ 
ration for a given energy E and 0{N^) operations for 
each product (i. e. G applied to several different vectors) 
provided the initial value of E is not changed. In terms 
of the matrix size N 2 « N'^/2 this implies a complex¬ 
ity of 0 {N 2 ^^) operations which is more expensive than 
the product H\(p> with 0 {N 2 ) operations but still much 
better than the naive matrix vector multiplication with 
G(iV|) operations. 

The position-energy transformation can be furthermore 
optimized for larger system sizes using that the one-particle 
eigenfunctions (j),y{x) are localized around some position 
x-max with localization length In this case the ratio 
\4'v{x)/(j),y{x^ax)\ is below 10“^^ (the numerical round¬ 
ing error for standard double precision numbers) for \x — 
a^maxl > c with the constant c = 17 log(10)£ « 175 if we 
replace the value ^ ~ 4.48 for A = 2.5. The positions x ful¬ 


filling this condition can be safely excluded in the multiple 
sums for the position-energy transformation therefore re¬ 
ducing the complexity to 0{cN^). 

This first variant of the algorithm combined with the 
(simple) Arnold! method for G is already very efficient and 
very superior to the implicitly restarted Arnold! method 
applied to H and produces for a sufficiently large value 
of the Arnold! dimension ua easily more than 50% — 70% 
of numerically accurate eigenvalues close to the energy E 
appearing in the Green function (from all ua Ritz eigen¬ 
values produced by the Arnold! method). For example for 
the Hubbard case with U = 7.8 and E = —2.78 we have 
been able, on a machine with 64 GB of RAM memory, 
to increase the system size up to A^ = 4181 (which is a 
Fibonacci number) and to choose the Arnold! dimension 
UA = 900 and about 620 out of 900 obtained eigenvalues 
have a quality with 6'^E{ip) < 10“^° [29|. Furthermore 
most of the important parts of the algorithm can be quite 
well parallelized for multiple core machines. 

As start vector for the Arnold! iteration we choose a 
vector proportional to the projection Pj^xi X 2 
i. e. a vector with uniform identical values for the sites in 
the set S where the interaction acts. In this way we avoid 
(most of) the many useless contributions from eigenstates 
which are essentially localized product states > 

with both particles localized very far away such that the 
interaction has no effect on them. With this start vector 
we capture all well “delocalized” states with energies close 
to the value of E. The Arnold! method still provides a con¬ 
siderable number of eigenstates being similar to strongly 
localized product states where the distance between par¬ 
ticles is “modest”, i. e. sufficiently large that the product 
states are indeed relatively good eigenstates of H but also 
sufficiently small that the initial vector has small contribu¬ 
tions of these states which will be amplified by the Green 
function Arnold! method if the eigenvalue of the product 
state is sufficiently close to E. 

For small values of Ur the memory requirement of the 
Arnoldi method is determined by the number ua of itera¬ 
tion vectors which need to be stored and the size of these 
vectors N 2 « N'^ /2 which provides the essential limita¬ 
tion of this method concerning the choice of ua and N. 
For larger values of Ur, e. g. Ur = 20 the largest value we 
have considered, the requirement to store multiple matri¬ 
ces of size UrNxUrN is also important (or even dominant 
for the second variant described below). 

However, this first variant of the Green function Ar¬ 
noldi method, which works with vectors stored in the posi¬ 
tion representation, can be considerable improved by using 
vectors stored in the non-interaction energy representation 
using an expansion in terms of the non-interacting prod¬ 
uct states \4>u,4>tJ. >■ This modification allows for several 
improvements. 

First, the number of the rather expensive energy-po¬ 
sition (or inverse position-energy) transformation steps is 
reduced from four to two when evaluating G\ip > since, 
according to the above description of the algorithm, the 
first energy-position and the last inverse position-energy 
transformation can be avoided if the vector |(/? > is by 
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default already available (or needed) in energy represen¬ 
tation (instead of position representation). 


Second, in this modified variant it is natural to choose 
a somewhat different start vector, i. e. a vector given as 
sum of product states with maximal positions in the set 
S which is qualitatively similar to the other initial vector 
used for the first variant but still different due to the finite 
one-particle localization length. The important point is 
that the new initial vector contains less contributions from 
useless products states. For given values of ua and N this 
improves considerably the quality of the eigenvectors by 
reducing the value of the quantity (26) and one obtains 
more nicely “delocalized” states (witn eigenvalues a bit 
further away from E) and less useless product states. 


The third improvement concerns the possibility to re¬ 
duce considerably the dimension of the Hilbert space in en¬ 
ergy representation from N 2 « N'^/2 to cN (with c « 175 
for A = 2.5) since one can simply remove all product states 
with maximal positions further away than c because these 
states do not feel the interaction at all (i. e. with interac¬ 
tion coupling matrix elements smaller than 10“^^). This 
reduces the amount of memory usage and also computa¬ 
tion time for the Arnoldi iterations by a factor 2c/N which 
becomes quite small for large system sizes {N > 1000). Es¬ 
pecially the reduced memory requirement allows to per¬ 
form computations with larger values of N and ua, for 
example for Ur = 1 we have been able to choose a system 
size N = 10946 with Arnoldi dimension ua = 3000 (on 
a machine with 64 GB of RAM memory). For the case 
N = 4181 and ua = 900, the maximum possible size for 
the first variant with 64 GB, the computation time for the 
second variant of the method is reduced by a factor of ten 
if compared to the first variant. 


The overall complexity of the Green function Arnoldi 
method for small systems {N < c) is given by Ci{UrN)^ + 
ua + C^N"^ n\ with three terms representing the ini¬ 
tial preparation part (first term with the constant Ci ~ 1), 
the Green function vector multiplications (second term 
with the constant C 2 5) and the Gram-Schmidt orthog- 
onalization scheme (third term with the constant Ca ~ 1). 
For larger systems N ^ c = 175 we have to replace in the 
second and third term a factor of A^ by c resulting in 
Ci{UrN)^ + C 2 cN'^ ua + CscN n\. If one choose typi¬ 
cally Ua N the second and third term have comparable 
complexity ^ c but in practice the second term is dom¬ 
inant due to a considerably larger value of the constant C 2 - 
Therefore it is not interesting to use the Lanczos method 
(instead of the full Arnoldi iterations) because this would 
only remove in the last, non-dominant, term one factor 
of N. The memory requirements (in units of size of dou¬ 
ble precision numbers) scale with C 4 {UrN)"^ -\- C^cN ua 
with C 4 ^ 5 C 5 and C 5 ^ 1 because one has to store sev¬ 
eral copies of matrices of size {UrN) x (UrN) and ua 
vectors of size cN for the Arnoldi iterations. 


With increasing values of the interaction range Ur the 
memory requirement and also computation time of the 
initial preparation part become more important or domi¬ 
nant, for Ur = 20, but even for this extreme case we have 
been able to push the system size up to = 1597 and 


one can (should) choose very large values for ua for the 
second Arnoldi-iteration part to better exploit the compu¬ 
tational “investment” of the preparation part. Even with 
UA = 2500 for N = 1597 the second and third part require 
only about 5% of the computation time while for Ur = 1 
the first preparation part is typically negligible (at most 
7% for the largest system size N = 10946, ua = 3000 we 
considered). 

We close this Appendix mentioning that the effective 
algorithm to compute arbitrary resolvent vector products 
can also be used to calculate more directly (or improve) 
individual eigenvectors if the eigenvalue (or an approx¬ 
imate eigenvector) is known with sufficient precision by 
the method of inverse vector iteration. We have for exam¬ 
ple been able to improve the modest quality eigenvectors 
which we had obtained by the implicitly restarted Arnoldi 
method to maximum possible precision only using a few 
number of these iterations. Actually, also a random initial 
vector can be used if a rather good approximate eigen¬ 
value is known. However, to achieve a good efficiency for 
a systematic computation of many eigenvectors with close 
energies the Arnoldi method for the resolvent is the best 
choice to exploit the Green function algorithm. The reason 
is the expensive initial part of the algorithm [the rather ex¬ 
pensive initial computation of Go and the matrix inverse 
in (|m] which is only done once for the Arnoldi method 
ananas to be repeated for any new individual eigenvalue 
when using inverse vector iteration. 


C Projected Green’s function formula 

In this appendix we show the formula (§ where G = {E — 
H) Gq = [E — Hq) H = Hq + U, Gq = PGqP 
and P = P^ is a projector such that tJ = PU = IJP = 
PUP, i. e. U has the same eigenvectors as P and only 
non-vanishing eigenvalues if the corresponding eigenvalue 
of P is unity. 


C.l Perturbative expansion of G 

The proof of (§ by an expansion in a matrix power series 
is quite illustrative. First we express G as 


G = 


il-UGo)iE-Ho) 


1 -1 


= Go(l-GGo)-i =Go^(GGo)^ 


= Go + Go ^(GGo)^ 


n—0 

’ ] UGo 


(31) 


(32) 


V,n=0 


where we have assumed that the matrix power series con¬ 
verges well which is the case for sufficiently large values 
of E in the complex plane. Using the relations between U 
and P we may rewrite the expression (32) as: 


G = Go + Go ^{UPGoPT UGo 


(33) 


\n—0 
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which becomes after replacing Gq = PGqP and resum¬ 
ming the series (in parentheses) just formula Further¬ 
more applying an argument of analytic continuation the 
validity of ®> is extended to all values of E in the complex 
plane (except the singularities of G or Go). This calcula¬ 
tion shows the crucial role of the relations between the in¬ 
teraction operator U and the projector P and which hnally 
allow to reduce the difficulty to determine the resolvent G 
by using a matrix inverse in a subspace of considerably 
smaller dimension which is just the subspace onto which 
P projects. 


C.2 Algebraic direct proof 


The expansion in a matrix power series and the argument 
of analytic continuation can be avoided by a direct but 
somewhat “less clear” calculation. For this we write: 


G = G{E - Hq)Go = G{E -H + U)Go = Gq + GUGq 
= Go + Go{l - UGo)-^UGo (34) 

= Go + GoOUGo (35) 


where we have used the first identity of ( |3l[ ) to obtain 
(34). The operator O is given by O = (1 — PA)~^P and 
A = UGq and to obtain (35) we have used (twice) that 
PU = U. We rewrite O in the form 


6 = {1-PA)-^P{1-PAP)il-PAP)-^ (36) 


and since P{1 — PAP) = (1 — PA)P we obtain the ex¬ 
pression 


6 = P(1 - PAP)-^ = (1 - PAP)-^P = (1 - UGo)-^P 


which together with (35) (and again PU = U) provides 
the formula (§. 
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